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PREFACE 


In February, 1965, Dr. Ernst Stuhlinger, now Marshall Space Flight 
Center's Associate Director for Science, initiated a series of Re- 
search Achievements Reviews which set forth those achievements 
accomplished by the laboratories of the Marshall Space Flight Center. 
Each review covered one or two fields of research in a form readily 
usable by specialists, systems engineers and program managers . 
The review of February 24, 1966, completed this series. Each re- 
viewwas documented in the "Research Achievements Review Series. " 

i 

In March, 1966, a second series of Research Achievements Reviews 
was initiated. This second series emphasized research areas of 
greatest concentration of effort, of most rapid progress, or of most 
pertinent interest and was published as "Research Achievements 
Review Reports, Volume II." Volume II covered the reviews from 
March, 1966, through February, 1968. 

This third series of Research Achievements Reviews was begun 
in March, 1968, and continues the concept introduced in the second 
series. Reviews of the third series are designated Volume III and 
will span the period from March, 1968, through March, 1970. 


Th* papmrm tn this r*porf warm prmamntmd Junm 26# 7969 


William G. Johnson 
Director 

Research Planning Office 


iii 




I 


CONTENTS. 

INTRODUCTION TO METEOROID PHYSICS RESEARCH 

By R. J. Naumann 

Page 

INTRODUCTION 1 

ACKNOWLEDGMENTS 3 

LIST OF TABLES 

Table Title ;• Page 

1. Energetics of Meteoroid Encounter 1 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Summary of various proposed meteoroid environmental models 2 

2. Current meteoroid environment 3 

^>1^HT GAS GUN EXPERIMENTS 

By David Jex 

Page 

FACILITY DESCRIPTION 5 

EXPERIMENTS CONDUCTED 5 

LIST OF TABLES 

Table Title Page 

1. Areas of Involvement and Specific Experiments Conducted . 6 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Simulated meteoroid impact on solar cell panel 6 

2. Simulated meteoroid impact of bumper shield configuration 7 


v 



CONTENTS (Continued). . . 


Page 


3. Bumper shield and secondary surface damage 8 

4. Results of bumper shield test series 8 

5. Comparison of normalized results of the bumper shield test series with a 

similar independent study 9 

6. Results of honeycomb material test series 10 

7. Relation of failure area to impact point for honeycomb material test series 10 

8. Results of calibration tests on Pegasus panels 11 

9. Output waveform profiles of Pegasus panels 11 

10. Target configuration for dielectric study 12 

11. Electromagnetic disturbance waveform associated with dielectric study 12 


-^SIMULATED METEOROID IMPACT TESTING ON MDA CONFIGURATION 


By R. C. Ruff 

Page 


SUMMARY 15 

INTRODUCTION 15 

HYPERVELOCITY LAUNCH TECHNIQUES 15 

PENETRATION TESTS 17 

CURRENT RESEARCH 25 

CONCLUSIONS 29 

REFERENCE 29 


LIST OF ILLUSTRATIONS 


Figure Title Page 

1. Materials Division Meteoroid Simulation Facility 16 

2. Schematic diagram of range tank 17 

3. Schematic diagram of improved sabot stripper 18 

4. Target configuration for MDA wall penetration tests 19 


vi 



CONTENTS (Continued). . . 

Page 

5. Impact test on MDA Series 1 19 

6. Basic MDA configuration compared both to magnesium bumper and to no insulation 20 

7. Basic MDA configuration compared to 0. 635-cm (0.250-in. ) backup 21 

8. Example of damage when no insulation was used 22 

9. Example of damage when Insulation was used 23 

10. Target configuration for MDA external electrical cable tunnel tests 24 

* 

11. Example of damage received by cable when protected only by bumper 25 

12. Target configuration for MDA window bumper tests .i 26 

13. Window condition immediately after test (left) and after being wiped with dry cloth (right) ... 26 

14. Window after cleaning with alcohol. . 27 

15. Multilayers of aluminum foil and high density foam 28 

16. Multilayers of aluminum foil and cork 28 

17. Multilayers of aluminum foil and red polyurethane foam 29 

18. Multilayers of aluminized mylar and red polyurethane foam 29 

19. Multilayer target density versus both penetrated thickness and areal density 30 

20. Bumper penetration as photographed by an image converter sequence camera 30 

21. Sequential views of a complete penetration of bumper and backup plate 31 

mO ACCELERATOR RESEARCH 

Page 

By Patrick N. Espy 33 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Plasma range showing firing console 34 

2. Plasma range from other side. 34 

3. 250 kJ capacitor bank 35 

4. Coaxial gun presently in use 35 


vii 



CONTENTS (Continued). . . 


Pago 


, 5. Plasma range with coaxial gun Installed 36 

6. Self-luminous plasma emerging from coaxial gun 36 


WlMPACT code development 

Page 

By R. J. Naumann 37 


LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Computer studies of impact phenomena represent the material by a series of imaginary 
mass points. As the impact progresses, the motion of these mass points indicates the 
material flow. The regions where the mass points cluster to form the darker areas are 
regions in which the material has been compressed by the strong shock waves generated 
by the impact. Such shock pressures can be as high as 100 million atmospheres (many 
times the pressure at the center of the earth). Material behavior at such pressures can 
be determined from quantum mechanical calculations and has been measured at more 
modest pressures (1 to 2 million atmospheres) by controlled shock experiments 37 


2. Dynamic pressure profiles 38 

3. Penetration versus impact velocity 39 


4. Decay of the maximum value of a at the axis with distance into the target 39 

zz 

5. Two-wave shock structure approaching pellet at rear surface of target 40 

6. Comparison of crater and splash for different size profiles at scaled times of 



Page 


INTRODUCTION 41 

THEORETICAL MODEL OF IMPACT-GENERATED SHOCK PROPAGATION 

IN SOLID MEDIA 42 

FORMULATION OF THE CHARACTERISTIC EQUATIONS . 43 

HUGONIOT RELATIONS OF SHOCK WAVES IN SOLID MEDIA 47 

SHOCK PROPAGATION LAWS 49 


vlii 



CONTENTS (Continued). . . 

Page 

SHOCK WAVE INTERACTION WITH FREE BOUNDARY 50 

CONSTRUCTION OF NUMERICAL SOLUTIONS FOR SHOCK WAVE PROPAGATIONS 51 

AN APPROXIMATION FOR PRESSURE PULSE AT LATE STAGE OF IMPACT 53 

NUMERICAL RESULTS AND DISCUSSION 55 

APPENDIX: SOUND SPEED OF c,s MATERIAL BEHIND THE SHOCK 59 

REFERENCES ‘ 60 

LIST OF TABLES 

Table Title Page 

1. Shock and Adiabatic Sound Speed Constants for Aluminum . ‘ 49 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. General features of the hypervelocity impact mechanism 42 

2. Moving shock front in a suddenly enlarged channel 42 

3. Ray tube coordinates system 43 

4. Geometric relation of equations (2) and (3) for the two-dimensional case 43 

5. Mach angle in physical plane 44 

6. Relation of (x,y) plane and (0 t t) plane 45 

7. Sketch of shock propagation in solid medium (t^ is the time when the first 
characteristic line reaches the shock centerline, x - axis; the normal shock 

portion vanishes completely for t= t^) 45 

8. Basic coordinate system of axisymmetric case 46 

9. Geometric relation for axisymmetric cases 46 

10. The ratio of adiabatic sound speed and the scope of Rankine-Hugoniot curve as 

a function of shock Mach number for c,s material 48 

11. Shock wave and free surface interaction 51 

12. Impact-generated normal shock waves 52 

ix 


CONTENTS (Continued). . . 


13. Construction of a simpio wave solution 

14. Approximate shock profile for the strong shock situation (U largo) 

15. Spherical coordinate system with symmetric axis 

16. The starting point of shock curvature in aluminum 

17. Numerical example of aluminum impact on aluminum at a velocity of 20 km/sec 

18. Comparison of peak pressure distribution for an aluminum-on-aluminum 

impact (impact velocity is 20 km/sec) 

19. Angular pressure distribution 


tfNE- DIMENSIONAL MODEL FOR INVESTIGATION OF CHANGES IN CONDUCTIVITY 
OF DIELECTRICS UNDER SHOCK CONDITIONS 


Page 

52 

53 

54 

56 

57 

58 
58 


By E. E. Klingman 

Page 

ABSTRACT 63 

INTRODUCTION 63 

MODEL FOR CALCULATIONS 63 

CONCLUSIONS 76 

REFERENCES 76 

LIST OF TABLES 

Table Title Page 

1. Clebsch- Jordan Coefficients 69 


LIST OF ILLUSTRATIONS 


Figure Title Page 

1. Formation of bands 64 

2. Energy bands for crystal with N atoms of ls 2 2s 1 configuration 64 

3. Four general types of band structure 65 


x 



CONTENTS (Continued). . . 

Page 

i 

4. An illustration of the difference in discrete and continuous energies of single 

atom and continuous energy bands and forbidden energies of periodic potential. . . 65 

5. Model chosen to determine Vj(o:) 66 

6. Coordinate system for ion-atom pairs 67 

7. Additive energy as a function of lattice parameter 68 

8. Slater representation of crystal potential 70 

9. Potential energy of electron in hydrogen and alkali atoms 70 

10. Variation of potential as a function of lattice parameter 71 

11. Illustrating n-dependence of bandwidth 72 

12. Energy band structure for model used 72 

13. Alternate method of depicting energy gaps resulting from model 73 

14. Different possibilities of overlapping wave functions 75 

TEMPERATURE EQUATION OF STATE FOR ALUMINUM 

By Ft. J. Naumann 

Page 

SUMMARY 79 

INTRODUCTION 79 

FORMULATION 80 

COMPUTATIONAL RESULTS 85 

CONCLUSIONS 92 

REFERENCES 93 

LIST OF TABLES 

Table Title Page 

1. Summary of Hugoniot Data for Aluminum 91 

xi 


CONTENTS (Continued). . . 


LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Comparison of solid and liquid potentials 83 

2. Behavior of heat capacity for solids and liquids in the high temperature regime 85 

3. Transition of heat capacity from liquid to gaseous behavior 86 

4. Gibbs energy isotherms 87 

5. Melting curve deduced from Gibbs energy plot 87 

6. Pressure isotherms 88 

7. Energy isotherms , 89 

8. Isotherms on S-V plot 89 

9. Isoenergy lilies 89 

10. D-u plots comparing computed results with experimental data points 90 

11. P-u plots comparing computed results with experimental data points 90 

12. T-u plots 92 

13. S-u plots 92 

LIGHT LEVEL TV OBSERVATIONS 

Page 

By Stuart Clifton 95 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Image orthicon and accessory equipment mounted in the DC -3 for airborne 

meteor observations 96 

2. A meteor as observed by the image orthicon system 97 

3. An angular rate distribution of meteors as recorded by the image orthicon 

during observations throughout November, 1969 97 

4. The decrease in apparent brightness of an object with increasing angular velocity 98 

5. A typical peak-height light curve as recorded by a photomultiplier tube observing 

the television monitor 98 

xii 



CONTENTS (Continued). . , 



PHOTOELECTRIC OBSERVATIONS OF METEORS TO 6.4 MAGNITUDE 


By U. J. Naumann 

SUMMARY 

INTRODUCTION 

THEORETICAL LIMIT OF A PHOTOMULTIPLIER METEOR DETECTION SYSTEM 

CALIBRATION 

METEOR OBSERVATIONS 

CONCLUSIONS 

REFERENCES 


Page 

101 

101 

101 

105 

106 
110 
no 


LIST OF TABLES 

Table Title Page 

1. Meteors Observed in 17 Hours, 22 Minutes Observing Time 109 


LIST OF ILLUSTRATIONS 

Figure Title Page 

1. 120-Hz sky background and level shift caused by a 3.45 magnitude meteor 

(50 mV/cm, 10 ms/cm) 107 

2. a. Typical meteors observed — 3.58-magnitude meteor (50 mV/cm, 100 ms/cm) 107 

b. 5. 16 -magnitude meteor (50 mV/cm, 50 ms/cm) 107 

c. 4. 5 -magnitude meteor (50 mV/cm, 50 ms/cm) 107 

3. Results of this study compared with photographic work of Hawkins and Upton 108 

4. Filter attenuation at different pulse durations . 108 

5. Influx rates of this study compared with those observed by Hawkins and Upton 110 


xlii 



CONTENTS (Concluded). . . 

OBSERVATORY PLANS AND ASTRONOMICAL APPLICATIONS 

Pago 

By James McGuire Ill 

LIST OF TABLES 

Table Title Page 

1. Observational Conditions of Observatory Site Ill 

LIST OF ILLUSTRATIONS 

Figure Title Page 

1. Location of observatory site at Redstone Arsenal Ill 

2. View of eastern horizon from observatory site 112 

3. View of southern horizon from observatory site 112 

4. View of western horizon from observatory site 113 

5. View of northern horizon from observatory site 113 

6. 12-foot diameter dome located at observatory site 114 

7. Ferson telescope having a 16 -inch primary mirror 114 

8. Celestron Pacific telescope 115 

9. Celestron telescope with a low light level TV camera and rich field adaptor attached 116 


xiv 



INTRODUCTION TO METEOROID PHYSICS RESEARCH 

By 

R. J. Naumann 


INTRODUCTION 

Since the time man first began to seriously 
consider spaceflight, it was recognized that the 
existence of small bits of material in space ranging 
from micron- size dust to huge boulders presented a 
possible hazard to space vehicles. This hazard 
results from the extremely high relative velocities 
of the particles, ranging from earth escape 
(11 km/sec) to parabolic retrograde (72 km/sec), 
with aa average of approximately 20 km/sec. The 
energetics of impacts at these velocities are sum- 
marized in Table 1. 

Many efforts have been made to assess the 
meteoroid hazard to space vehicles. The problem 
is two-fold. First, the natural environment must be 
defined in terms of flux or number density as a 
function of mass and velocity. Secondly, the damage 
to a given structure as a function of mass and velocity 
must be specified. With these two ingredients, the 
hazard to a particular space vehicle may be deter- 
mined by integrating over the mass and velocity 
distributions to find the probability of encounter with 
a particle that will cause critical damage. 


Early attempts to define the mass -flux distri- 
bution relied on visual observation of meteors and 
subjective estimates of their brightness compared to 
fixed stars. Later, photography provided a more 
precise method for obtaining magnitudes. However, 
photographic methods are limited to about first or 
second magnitudes. The problem of interpreting the 
observed brightness in terms of energy or mass 
remained. This was solved by assuming some 
luminous efficiency, or fraction of energy converted 
to light. Later experiments such as Trailblazer, 
which launched an artificial meteor of known mass 
and velocity into the atmosphere, provided a some- 
what firmer basis for this interpretation. 

Structural damage was estimated by extrapolating 
ballistic data at a few kilometers per second to 
meteoroid velocities — approximately a five -fold 
extrapolation. 

The number of hazard models was almost equal 
to the number of investigators, since each had his 
own opinion of how to do the various interpretations 
and extrapolations, as may be seen in Figure 1. It 
was not until penetration experiments were flown 


TABLE i. ENERGETICS OF METEOROID ENCOUNTER 



1 1 km/ sec 

20 km/ sec 

72 km/ sec 

Available Energy (MJ/kg) 

60 

200 

2500 

(eV/ amu) 

0.6 

2 

25 

Typical Pressures N/m 2 

50 x 10 9 to 30 x 10 10 
(0.5 to 3.0Mb) 

15 x 10 10 to 75 x 10 10 
(1.5 to 7.5Mb) 

15 x 10 11 to 75 x 10 11 
( 15 to 75 Mb) 

Temperature of Com- 
pressed Material (°K) 

0.5 x 10 5 

1.5 x 10 s 

2 x 10 6 


FOR COMPARISON: 


Yield of TNT is 1. 75 MJ/kg 

Typical metallic bond is 0.073 eV/amu (Fe) 


Largest ionic bonding energy is 0.4 eV/amu (LiF) 
Pressure at interior of earth is 4 x 10 11 N/m 2 (4.0 Mb) 
Temperature at interior of sun is 15 x 10 € *K 
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Figure 1. Summary of various proposed meteoroid environmental models. 

on satellites that some coherent picture began to Considerable improvements have been made in 

emerge. Through use of the more recent measure- the development of hypervelocity techniques as well 

ments and a more refined interpretation, the present as in theoretical understanding of impact phenomena 

picture of the meteoroid environment is illustrated and of the behavior of matter at the pressures and 

in Figure 2. temperatures associated with meteoroid impact. 
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( PARTICLES PER SQUARE METER-SECOND OF MASS, m, OR GREATER ) 


ft. J. NAUM ANN 



MSFC has contributed heavily in the evolution of 
the present picture of the meteoroid environment as 
well as in experimental and theoretical work in 
accelerator development, impact mechanics and 
related phenomena, and high-pressure, solid-state 


physics. It is the purpose of this review to report 
on some of the more recent accomplishments. 

ACKNOWLEDGMENTS 


The Astrophysics Branch of the Physics & 
Astrophysics Division, Space Sciences Laboratory 
(SSL) wishes to express its gratitude to various 
NASA Headquarters agencies, expecialiy to OART for 
its support of the meteoroid program at Marshall 
Space Flight Center (MSFC) . We particularly wish 
to thank Mr. Charles D’Aiutolo for his interest and 
encouragement in this program . Also we are 
deeply indebted to SSL personnel, particularly our 
Director, Mr. Gerhard B. Heller, and to our former 
Director, Dr. Ernst Stuhlinger, and to our former 
Division Chief, Dr. James B. Dozier, and his deputy, 
Mr. Ray V. Hembree, for their strong endorsement 
of this activity and for providing management support 
to these programs . 

We wish to express our appreciation to the other 
MSFC Laboratories that supported this effort and in 
particular to Mr. Clyde Nevins, of S&E-ASTN-ES, 
who has cooperated with us extensively in developing 
design information. 

We also wish to acknowledge the efforts of our 
technicians, Mr. Archie Miller, Mr. Huie James, 
and Mr. Sam Brown of Brown Engineering Company, 
and the efforts of Mr. Charles MacKay and Mr. Jim 
Harrison of S&E-SSL-P, our summer student, William 
Peters, and our co-op students, Mike Duggan, David 
Rannefeld, Jerry Washburn, William Langley, 

Jennifer Hornbuckle, Wanda Neely, and Dean Mitchell, 
all of whom spent many hours observing and analyzing 
data. 


3 





k N?0-35 oO 6 

LIGHT GAS GUN EXPERIMENTS 

By 

David Jox 

FACILITY DESCRIPTION meteoroid detectors and basic research for future 

applications in meteoroid detection analysis. 


At present the Space Sciences Laboratory (SSL) 
is operating two light gas guns. The smaller gun 
has a nominal 0. 1588-cm (0.0625-in.) launch tube, 
but is capable of launching projectiles up to 0.2381 cm 
(0.0938 in.) in diameter. Driven by a .375 H&H 
Magnum rifle cartridge, the gun has a velocity 
capability of 8.5 km/sec. The larger gun has a 
nominal launch tube diameter of 0.318 cm (0. 125 in.) 
that can be extended to 0.635 cm ( 0. 25 in. ) . A 
20 -mm aircraft cannon cartridge is used as a driver. 
Velocities approaching 10 km/sec have been achieved 
A novel sabot ing technique allows glass beads as 
small as 100 \x as well as full bore projectiles to be 
launched. 

Instrumentation for these ranges includes 
standard instrumentation such as image converter 
cameras, 100 -MHz counters, photomultiplier tube 
velocity stations, etc. Also there are two field 
emission Febetron 600 keV pulsed electron guns 
available that can be used to produce betagraphs, 
flash x-rays, or super -radiant light. 

The primary intended uses of this range were to 
study physical phenomena associated with hyperveloc- 
ity impact and to develop and calibrate meteoroid 
detectors. However, the high velocity capability and 
low operating cost of our range made it very much in 
demand for providing quick answers to actual design 
problems at MSFC and at other NASA installations. 

EXPERIMENTS CONDUCTED 

The light gas gun experiments conducted in the 
last two years can be grouped under three main 
headings. 

1. Assistance in quick analysis of basic struc- 
tures exposed to the meteoroid environment. 

2. Examination and calibration of existing 


3. Basic research in understanding the process 
and characteristics of the impact phenomenon. 

These three main headings together with those 
experiments that logically fall under each are shown 
in Table 1. The eight specific tests listed under 
these headings will be discussed in some detail in 
the following pages. 

Assistance in Quick Analysis of Basic 
Structures Exposed to the Meteoroid 
Environment 


VULNERABILITY OF SOLAR CELLS 

At the request of the Astrionics Laboratory 
of MSFC, the vulnerability of solar cells to be used 
on the Apollo Telescope Mount was examined. The 
questions to be answered are: Will a panel of solar 
cells suffer a catastrophic failure when impacted by 
a meteoroid and what size meteoroid can the panel 
withstand and still not suffer catastrophic failure? 

Figure 1 shows the results of an impact of a 

0. 3 18 cm (0. 125 in, ) diameter, 0. 318 cm (0 . 125 in. ) 
long cylindrical projectile with a density of 

1. 2 gm/cm 3 at a velocity of 7. 92 km/sec. By 
visually examining the extent of damage and measuring 
the voltage produced under the same conditions by 

the panel before and after impact, it was evident that 
only those cells directly impacted failed. There was 
no evidence of catastrophic failure such as propa- 
gations of cracks or delaminations. Since the pro- 
jectile size used in this test represented an upper 
limit on what could be expected on a large solar cell 
array, it was concluded that large meteoroids did 
not present a severe hazard to solar cell arrays. 
Additional tests with smaller projectiles and at 
oblique impacts defined the expected failure rate of 
solar cell panels. 
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TABLE 1. AREAS OF INVOLVEMENT AND SPECIFIC EXPERIMENTS CONDUCTED 


1. ASSISTANCE IN QUICK ANALYSIS OF BASIC STRUCTURES 
EXPOSED TO THE METEOROID ENVIRONMENT 

a. Vulnerability of Solar Colls 

b. Bumper Shiold Effectiveness 

- c. Protective Characteristics of Honeycomb Material. 

d. Evaluation of Meteoroid Shield on Mariner 71 
(Jet Propulsion Laboratory) 

e. Check on Soft X-ray Device (Naval Research Laboratory) 

2. EXAMINATION AND CALIBRATION OF EXISTING METEOROID 
DETECTORS AND BASIC RESEARCH FOR FUTURE APPLICATION 
IN METEOROID DETECTION ANALYSIS 

a. Calibration of Pegasus Panels 

b. Electromagnetic Pulse Associated With Shocked DieLectrics 

3. BASIC RESEARCH IN UNDERSTANDING THE PROCESS AND 
. CHARACTERISTICS OF THE IMPACT PHENOMENON 

a. Research of Impact Phenomenon (Contract With 
IIT Research Institute) 



Figure 1. Simulated meteoroid impact on solar cell panel. 
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BUMPER SHIELD EFFECTIVENESS 

When it was determined that additional meteoroid 
protection was needed for the Orbital Workshop, it 
was proposed to deploy a thin shield or meteoroid 
bumper around the cylindrical S-IVB tank. To deter- 
mine the effectiveness of this bumper, it was nec- 
essary to determine what size projectile could per- 
forate this shield and remain relatively intact. 
Empirical data to answer this question were generated 
in a test series in cooperation with the Astronautics 
Laboratory, MSFC. The approach taken by this lab 
was to fix the spacing between the bumper or shield 
and the secondary surface representing the actual 
spacecraft surface. Then, by varying the projectile 
size and bumper thickness, find some relations that 
could aid in extrapolating to the expected meteoroid 
impact condition. 

Figure 2 shows the actual impact process as 
viewed by the TRW image converter camera. The 
profile of the debris cloud can be seen. Figure 3 
is a photograph of one of the specimens used. The 
damage pattern is typical of the results using this 
cylindrical type projectile. This particular pro- 
jectile impacted at an angle as may be seen by the 
oblong hole in the bumper. This also produced the pec- 
uliar spray pattern on the rear surface. Projectiles 


impacting end-on produced the conventional 
circular spray pattern. In either case, the depth 
of penetration into the rear surface is essentially 
the same. Figure 4 is a graph of the results. It will 
be noted that the depth of damage in the secondary 
surface remains constant for each size projectile 
until it reaches a point where the bumper is no longer 
sufficient to completely break up the projectile. At 
this point, there is a noticeable increase in the depth 
of the damage. The minimum bumper is the thinnest 
shield necessary to completely fragment the projectile. 
Any thicker bumper will give adequate protection as 
far as depth of penetration is concerned. Figure 5 
is a comparison of the normalized data that were 
obtained with lexan on aluminum compared with 
similar work with aluminum on aluminum by C. J. 
Maiden at General Motors Defense Research Labo- 
ratory. From the figure, It is evident that the bumper 
is more effective against the lower density lexan 
(p = 1.2 gm/cm 3 ) than against aluminum with 
p = 2. 77 gm/cm 3 . This effect could be very significant 
since the bulk density of meteoroids is believed 
to be less than 1.0 gm/cm 3 . 

Additional information on this test series was 
presented at the OART Meteoroid Impact and 
Penetration Workshop October 8-9, 1968, at the 
Manned Spacecraft Center. 



Figure 2. Simulated meteoroid impact of bumper shield configuration. 
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Figure 5. Comparison of normalized results of the bumper shield test series with a similar 

independent study. 


PROTECTIVE CHARACTERISTICS OF HONEYCOMB 
MATERIAL 

Honeycomb material offers attractive advantages 
as a structural material because of its high rigidity - 
to-weight ratio. Also, the outer wall naturally serves 
as a meteoroid bumper, which should improve its 
protection capability over an equivalent weight single 
wall. The question arose as to whether the channeling 
action of the honeycomb structure might offset some 
of the gain afforded by the double wall. A series of 
tests was conducted to evaluate this question. 

The tests were conducted by holding projectile 
mass constant and varying velocity until failure of the 
rear surface was achieved. Figure 6 illustrates the 
results of these tests. The cutoff point near 80 
degrees was estimated from impacts at 75 degrees 
to the normal. Figure 7 shows the area of failure 
in relation to the point and angle of impact. It will 
be noted here that the failure area is within the 
perimeters of the cell of initial impact, although as 


indicated by the dashed lines, the cell walls were 
destroyed in most cases. This shows that the 
honeycomb directs, or channels, most of the debris 
resulting from the primary impact at a point. 

Tests on equivalent double-wall structures have 
not been completed at this time and will be published 
in a forthcoming report. 

EVALUATION OF METEOROID SHIELD ON MARINER 
71 AND CHECK ON SOFT X-RAY DEVICES 

Because of the capability to launch very small 
glass beads at high velocities, MSFC was asked by 
several outside agencies to assist in meteoroid 
testing of space vehicles and experiments. A series 
of shots was conducted at the request of the Jet 
Propulsion Laboratory to evaluate the meteoroid 
shield effectiveness of the insulation on Mariner 71. 
Tests were also conducted on a soft x-ray experi- 
ment under development by the Naval Research 
Laboratory. 
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Figure 6. Results of honeycomb material test series. 



Figure 7. Relation of failure area to impact point 
for honeycomb material test series. 

Examination and Calibration of Existing 
Meteoroid Detectors and Basic Research for 
Future Applications in Meteoroid 
Detection Analysis 

CALIBRATION OF PEGASUS PANELS 

The Pegasus spacecraft carried the first 
large area meteoroid detectors into space. The area 


of surveillance afforded by these spacecraft allowed 
a sampling of the meteoroid population in the size 
range from which significant structural damage can 
be expected. These data points form the basis for 
extrapolation of the meteoroid environment to the 
larger masses that are of critical interest to large 
manned vehicles. Since these are the only data 
concerning thiB size range that are expected to be 
available in the forseeable future, a substantial 
effort was made to calibrate the detector. 

Figure 8 shows the results for the 0. 04-cm 
Pegasus panel. Each vertical line represents a 
mass. The hash marks on the right side of the line 
represent detected impacts, while the hash marks 
on the left side of each line represent undetected 
impacts. In this manner, the threshold velocity 
may be established for a number of masses. 
Extrapolation of these points to average meteoroid 
velocity gives the average threshold mass for the 
detector. It may be seen that the higher velocities 
attainable by the SSL facility extended the calibration 
data to substantially smaller masses and thereby 
increased the confidence in the extrapolation. Also 
it is interesting to compare these results with the 
predicted results using the empirical penetration 
formula developed by Richard Fish at Ames Research 
Center. It may be seen that there is even a greater 
velocity dependence than that predicted by Fish. 

Figure 9 shows some of the output waveforms 
from the Pegasus detector. The photographs are 
those taken of the oscilloscope that monitored the 
panel discharge. The photographs illustrate 
the difference in penetration and nonpenetration. The 
typical sharp quick discharges on the left are those 
associated with penetration. The relatively long 
discharges or no discharges displayed on the right 
are associated with nonpenetration firings. The 
interesting fact is that a sharp deep discharge is 
obtained when the impact is substantially above 
threshold. However, when near threshold conditions 
prevail, the discharge becomes more gradual and 
more shallow until it is no longer detectable. 

ELECTROMAGNETIC PULSE ASSOCIATED WITH 
SHOCKED DIELECTRICS 

Several interesting phenomena associated with 
impacts on dielectrics have been observed recently 
in the form of electromagnetic disturbances. An 
insulated copper plate was installed behind a 
dielectric sheet. A shielded box surrounding the 
plate served as ground reference as shown in Figure 
10. A nonpenetrating impact on the dielectric 
produces a short positive pulse followed by a long 
negative deflection, as shown in Figure li. The top 
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Figure 9. Output waveform profiles of Pegasus panels. 
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Figure 10. Target cohfigu ration for dielectric study. 

profile is at 10 V/cm and the bottom profile is at 
0. 05 V/cm. The time base is the same for both. 
This long term negative deflection is obviously a 
result of a charge on the dielectric sheet that may 


have been deposited by the projectile, or may have 
resulted from chargo ojectlon during the impact. 

The origin of the fast |x>sllivo spiko is not so apparent. 
Additional tesLs will bo |>orformod to investigate thiH 
phenomenon and to determine Its potential ns a 
meteoroid impact detector. 


Basic Research in Understanding the 
Process and Characteristics of the 
Impact Phenomena 

RESEARCH ON IMPACT PHENOMENA 

To make certain that the basic phenomena 
associated with hypervelocity impact are understood, 
an experimental program has been conducted to make 
basic measurements on the debris cloud that can be 
compared with computed values. This program has 
concentrated primarily on the penetration products 
emerging from a thin sheet when impacted by a flat 
disk. This configuration was chosen to minimize 
jetting effects. Also the penetration products at the 
rear of a thin shield are of considerable importance in 
analyzing bumper performance. 



Figure 11. Electromagnetic disturbance waveform associated with dielectric study. 
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Most of the experimental work was carried out 
under contract NAS8-20337 at the Illinois Institute of 
Technology Research Institute by Jerome Jeslis. 
Measurements include velocity and density distri- 
butions in the debris cloud, temperature estimates 
based on photometry and spect ral data, time 


resolved spectra, and ionization state. All measure- 
ments were done at less than 13.3 x 10 N/m 2 ( 10 
torr) to minimize the effect of collision with ambient 
atmosphere molecules. The measurements are 
almost completed and will be published in a contractor 
report. 


13 




t N70-35oO 7 

SIMULATED METEOROID IMPACT TESTING 
ON MDA CONFIGURATION 


R. C. 


SUMMARY 

The results of research programs to improve 
the experimental capabilities of the Materials Divi- 
sion Meteoroid Simulation Facility of the Astronautics 
Laboratory are presented herein. The results of 
penetration threshold measurements on various 
configurations of the Multiple Docking Adapter (MDA) 
such as the spacecraft wall, external electrical cable 
tunnel, and photographic window bumper are dis- 
cussed. Also given are results of penetration tests 
on various densities of multilayer insulation mate- 
rials. The current research efforts are described. 

INTRODUCTION 


The programs that are presented in this paper 
were carried out in the Materials Division Meteoroid 
Simulation Facility of the Astronautics Laboratory. 
The programs had a dual purpose. One was to give 
basic data for design purposes and structural analysis. 
The second was to arrange the tests so that a good 
understanding of the physics of hypervelocity impact 
could be acquired. 

The following topics will be presented; ( 1) a 
brief summary will be presented of the research on 
hypervelocity launch techniques that were necessary 
to provide the experimental capabilities needed for 
penetration testing, (2) some of the results of the 
MDA configuration tests will be presented, (3) the 
results of a series of low density multilayer insulation 
tests will be related to protection efficiency , and 
(4) some current research studies will be briefly 
described. 

HYPERVELOCITY LAUNCH TECHNIQUES 

The Materials Division Meteoroid Simulation 
Facility is shown in Figure 1. The facility consists 


Ruff 

of a two-stage light gas gun with a 1 . 27 cm ( 0 . 5 in . ) 
diameter launch tube. The breech of the light gas gun 
is on the left with the pump tube, high pressure 
section, and launch tube extending to the right. At 
the far right is the dump tank for dispersing the 
acceleration gases, the x-ray diagnostics section for 
velocity and integrity determinations, and the target 
tank. 

This light gas gun facility was developed for the 
Materials Division under a research contract with 
the intent of providing an inherent flexibility in 
testing procedures. The success of this endeavor is 
the fact that projectiles with a mass range from 
0.008 to 5.0 grams can be launched at velocities 
from 5 to 8 km/sec. This facility was delivered 
approximately 2 years ago. After the first operational 
checkout firings, it was felt that there were two areas 
associated with the launch techniques in which 
improvements would greatly increase the information 
output of the facility. 

The first problem area was in the design of the 
sabot stripper that was being used. The sabot 
technique is a method of launching projectiles of 
various sizes and shapes. The projectile of interest 
is placed on the front face of a full bore, polycar- 
bonate slug called the sabot. The sabot and projectile 
must then be separated after leaving the muzzle of 
the light gas gun. This separating or stripping 
operation can be described using the schematic 
diagram of the range tank shown in Figure 2. The 
projectile and sabot (not shown) are launched from 
the barrel of the light gas gun shown on the left. As 
they leave the barrel, the sabot impacts four tungsten 
stripper pins, while the projectile, being less than 
full bore, goes between the pins and straight down 
range through a hole in the sabot trap and onto the 
target. The impact between the sabot and the stripper 
pins is a hypervelocity impact and, therefore, 
completely demolishes both the sabot and the pins. 

The debris from this impact is then caught on the 
sabot trap. 
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Figure i. Materials Division Meteoroid Simulation Facility. 


It was noticed that with the original design of 
the stripper pins some debris was getting through the 
hole in the sabot trap and also impacting the target. 
Although in most cases the debris damage could be 
separated from the projectile damage, there were 
several cases where the results were ambiguous 
and the test had to be repeated. To avoid this 
repetition, a short research effort was launched to 
try to understand the physics of shock-wave focusing 
and momentum transfer in this type of hypervelocity 
impact. The parameters available for studying 
changes in stripper design were the debris pattern 
on the sabot trap and the physical condition of the 
projectile. 

This research effort resulted in a new stripper 
design shown in Figure, 3. It was found that the front 
face of the tungsten piris had to be machined with 


approximately a 15-degree angle. Too large an 
angle focussed the shock wave to such a degree that 
the projectile was caught and damaged. Too small 
an angle meant that there was not sufficient radial 
momentum transfer, and some of the dense tungsten 
debris would go straight down range through the trap 
and impact the target. It was also found that the 
leading face of the sabot had to have a bevel that 
matched that of the stripper pins. If the face of the 
sabot is left flat, then the impact angle between the 
pins and the sabots creates a high velocity jet of 
material that is forced inward and is sufficient to 
damage the projectile or actually force it off axis so 
that it will hit the edge of the sabot trap. 

The evidence of the success of this new design 
is that in the last 120 shots there has not been a 
single failure because of the stripper operation. A 
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negligible amount of debris gets through the trap 
and the projectile always impacts within 
1 . 27 cm ( 0 . 5 in. ) of the center of the target. 

The second problem area was in velocity 
prediction. It was found that the velocity of the 
projectile could not be controlled to an accuracy of 
better than 8 or 10 percent. This is a common prob- 
lem in many facilities. The standard solution is to 
take a statistical distribution of velocities to cover 
the range of interest. It was felt that there could be 
a considerable savings in money and especially time 
in making penetration threshold determinations if the 
velocity could be predicted more accurately. 

For these reasons, a program was instituted to 
measure the variables in the first stage of the light 
gas gun such as the breech pressure, propellant 
burning time, and piston velocity. After analyzing 
these measurements, several changes in the loading 
and firing procedures were instituted. The velocity 
uncertainty is now down to ± 1 percent, which is 
essentially the accuracy of our velocity measuring 
instru mentation. 

After these two design improvements were insti- 
tuted, it was felt that the facility was ready to begin 
penetration threshold testing. The facility was ready 
to launch projectiles of controlled mass and velocity 
with negligible debris or velocity uncertainty. This 


meant that the penetration threshold for an unknown 
configuration could be zeroed in on rapidly and 
accurately. 

PENETRATION TESTS 
MDA Wall Configuration 


The MDA wall configuration was tested to 
determine its penetration threshold. The basic 
configuration is shown in Figure 4. This specific 
configuration consists of a 0.0508 cm (0.020 in.) 
thick 2014-T6 aluminum bumper on the left, a 
5.08-cm (2.0-in.) void, a 2.54 cm (1.0 in.) thick 
section of high performance multilayer insulation, and 
a 0. 196 cm (0.077 in.) thick 2014-T6 aluminum back- 
up plate that simulates the MDA pressure wall. The 
multilayer insulation consists of many layers of 0. 25 
mil doubly aluminized Mylar with red polyurethane 
foam spacer material. 

This basic configuration was tested to determine 
its penetration threshold as a function of mass and 
velocity. The results are shown in Figure 5 where 
the log of the mass is plotted versus the log of the 
velocity. The test points where no penetration occurred 
are shown in the filled circles, while penetration 
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Figure 3. Schematic diagram of improved sabot stripper. 


points are shown with open circles. The solid line 
shows the maximum slope that the penetration thresh- 
old curve could have over this velocity range. This 
is compared to the dashed line which is the Fish curve 
[ 1] , a penetration curve for a single sheet of alumi- 
num. 

Comparison between the two curves shows that 
the experimental configuration is much more efficient 
at defeating higher velocity projectiles than is a single 
sheet of aluminum. Another way of comparison would 
be to extrapolate both curves to meteoroid velocities. 
One then would find that the experimental curve gives 
a much better safety factor than the Fish curve. 


Unfortunately, without an experimental method of 
reaching at least 12 km/ sec, the extrapolation of the 
experimental curve is presently unjustified. There- 
fore, structural design engineers must use the Fish 
curve extrapolation to calculate the penetration 
resistance at meteoroid velocities. Using this 
calculation, it turns out that the basic configuration 
is equivalent to approximately 1.27 cm (0.5 in.) of 
single sheet aluminum. 

Penetration threshold determinations were made 
on several configurations that were variations of the 
basic configuration. The results from two variations 
are shown in Figure 6. The first variation was to 
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Figure 4. Target configuration for MDA wall 
penetration tests. 



Figure 5. Impact test on MDA Series 1. 


replace the 0.0508 -cm (0.020 -in.) aluminum bumper 
with an equal mass 0.08i3-cm (0.032-in.) magnesium 
bumper. These results are the open penetration 
triangle and the closed, no -penetration triangles. 

These points again delineate a line of maximum 
slope. This is the dashed Line slightly above the 
solid line which is the curve for the basic penetration 
from Figure 5. Therefore, it is seen that the 
magnesium bumper is more effective in defeating 
a hypervelocity projectile than an equal mass but 
more dense aluminum bumper. 

.The second variation shown in Figure 6 is the 
removal of the insulation from the basic configuration. 
Another dashed line showing the maximum slope is 
drawn through the square data points. It can be seen 
that without insulation the configuration is penetrated 
by a projectile of half the mass of that which is needed 
to penetrate the basic configuration with insulation. 
This curve, therefore, shows the great effectiveness 
an external multilayer insulation has in increasing 
penetration resistance. Also, it should be noted 
that the curve without insulation is essentially parallel 
to that with insulation. 

Test results on other variations of the basic 
MDA configuration are shown in Figure 7. These 
tests were made on targets with 0.635 cm (0.250 in. ) 
thick backup plates instead of the basic 0. 196-cm 
(0.077-in. ) plate. This configuration is also planned 
for use on specific parts of the MDA. These tests 
were also made with insulation (circles) and without 
insulation (squares). Tests on the targets with insu- 
lation provided a penetration threshold point. How- 
ever, the four data points were not sufficient to pro- 
vide a threshold curve of maximum slope. The data 
point for targets without insulation did give a 
maximum slope. 


The results shown in Figure 7 can be analyzed by 
comparison to the solid line, which is the curve for the 
basic configuration from Figure 5, and to the Fish 
curve for single sheet. Indications are that for such 
thick backup plates the threshold curves are more 
nearly parallel to the Fish curve. This is reasonable 
when it is noted that a 0.30-gram projectile is 
approximately 0.635 cm (0.250 in.) in diameter. 
Therefore, the bumper thickness to projectile diameter 
ratio is less than 1: 10. This is much less than the 
optimum range for efficient projectile breakup. This 
means that this configuration may act very much like 
a single sheet since the bumper does not contribute 
much protection. One further observation is that the 
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Figure 6. Basic MDA configuration compared both to magnesium bumper and 

to no insulation. 


addition of insulation to the configuration does not 
double the mass protection as it did at lower masses. 

As further evidence of the effectiveness of the 
multilayer insulation, it is interesting to compare 
the damage mechanism in a target without the insu- 
lation to one with insulation. Figure 8 shows the 
damage on a backup plate without insulation. The 
damage is a result of discrete particles coming off 
the bumper. Each particle is capable of cratering 
the backup plate, and where there is a sufficient 
density of particles, they actually knock holes through 
the backup plate. 


Figure 9 shows a target in which insulation was 
used. Approximately the same projectile and bumper 
were used as in Figure 8; however, total damage 
comparisons should not be made because different 
backup plate thicknesses were used. In Figure 9 it 
can be seen that there is no longer any particulate 
type damage. What has happened is that the insulation 
completely vaporized the small particles that were 
formed by the bumper impact. Therefore, the damage 
with insulation is a result of a rapid, high pressure 
vapor pulse that deformed the plate beyond its elastic 
limit and actually ruptured the backup plate along 
radial lines. 
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MDA Electrical Cable Tunnel 

Penetration tests were also made on the MDA 
external electrical cable tunnel. As the name implies, 
external cable on the MDA may also need to be 
protected against meteoroid damage. Several possible 
protection configurations were tested, but only the 
damage mechanism will be presented here. 

Figure 10 shows the test fixture used for these 
tests. The target consists of a 0. 102 -cm(0. 040-in. ) 
aluminum bumper that shields four bundles of 
electrical cable. Figure 11 shows the damage that 
is done to the cable during a simulated meteoroid 


impact. It can be seen that the debris from the 
bumper was in particulate form, as evidenced by the 
specific cratering on the backup plate. This partic- 
ulate matter also impacted on the electrical cable 
and vaporized much of the electrical insulation. How- 
ever, there was very little actual damage to the cables 
themselves. Also, after the impact it was seen that 
none of the stripped cables was physically touching. 

A possible reason for this is that when the debris 
impacted the wires, the gas pressure generated by 
the vaporization of the electrical insulations was 
sufficient to force the cables apart. 

Although there is no physical contact between 
cables, there are still transient effects that need to 
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Figure 8. Example of damage when no insulation was used. 


be considered. The energy in any hypervelocity 
impact, such as that between the bumper debris 
and the electrical insulation, is of sufficient magnitude 
to form an electrically conductive plasma. Therefore, 
during an impact, even without physical contact, 
electrical arcs will be possible between any 
electrical cables with a voltage difference. For this 
reason, several double-bumper configurations were 
tested but will not be discussed here. 


MDA Photographic Window Bumper 

Penetration tests were also made on the MDA 
photographic window bumper configuration. This 
configuration is shown in Figure 12. This config- 
uration consists of a 0.318 cm (0. 125 in.) thick 
fiberglass laminate bumper with 2.54 cm (1.0 in. ) 
of the high performance multilayer insulation bonded 
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Example of damage when insulation was used. 


to the bumper. In this configuration the 

5 08-cm (2.0-in.) void is between the insulation 
and the ultraviolet grade glass of which only a small 
is visible. The position of insulation is 
obviously for the purpose of swinging both the bumper 
and insulation aside when photographic experiments 
are desired. 

The penetration tests were made with two 
goals in mind. First, it was desired to f^ ^e pene- 
tratlon threshold for the bumper and insulation o y. 
Secondly, it was desired to determine ^tdamage 

is done to the window when penetration o the bum^r 

and insulation does occur. The penetration tests were 
made at approximately 6.8 km/sec . To find the 
threshold, the projectile mass was varied from 0.045 
<rr am down to 0.008 gram. Although the vapor 
pressure pulse from the impact of the 0 . 008 -gram 
projectile on the bumper was sufficient to tear a 
small slit through all layers of 
was considered the approximate threshold mass 


because no particulate material reached the window 
to cause damage. 

Figure 13 shows the damage that was done to 
the 1.2? cm (0.5 in.) thick window g^ssby a 

0 018 -gram projectile at 6.5 km/ sec. 

on the Eft shows the window before it was cleaned 
in any way. It can be seen that the window is very 
dar/ It is covered with the fine powder Produced 
hv the pulverization of multilayer insulation by the 

debris. This powder adheres quit, well a 

is evidenced by the view on the right wh.^hows the 
same window after it has been ^wtped I off a Uy 

paper towel. Some of the surface chippie* : »• 
caused by the particulate debris impact can be 

Figure 13 also, it ehouid be e^.a^d «« 
even with the larger 0 . 045 -gram projecUle ^ 
chipping was strictly surface damage and no structur 

damage was evident. 

Figure 14 is a close-up view of the same glass 
target Sown in Figure 13. This view was taken 
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Figure 10. Target configuration for MDA external electrical cable tunnel tests. 


after the glass had been cleaned with alcohol. It 
can be seen that all of the black powder has been 
removed. There are no smudges on the polished 
front surface to obscure the saw cut grain on the 
back surface. The small chips on the glass surface 
are also quite visible. Although hard to see in this 


view, the largest chip shows the conch oidal chip- 
out, which is the normal damage mechanism in a 
hypervelocity impact on glass. 

The results of these tests show that the window 
can be protected from meteoroid impact damage 
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Figure li. Example of damage received by cable 
when protected only by bumper. 

from a structural standpoint. However, meteoroid 
damage of the magnitude tested will roughen the 
surface by chipping over a localized area so that the 
camera would have to be displaced to an undamaged 
part of the window. Even this would not help until 
the black pulverized powder is cleaned from the out- 
side of the window with a solvent. 


Low Density Multilayer Insulation 


Penetrations were also made on a series of 
multilayer insulation targets of varying density to 


determine what effect density has on the efficiency 
of defeating a projectile. These targets were all 
tested with a 0.0 10 -gram aluminum projectile 
traveling at 6. 1 km/sec. 

Figure 15 shows the most dense target that 
was tested. It consists of aluminum foil sheets 
spaced with a relatively high density adhesive -backed 
foam. The overall density was 512 kg/m 3 . It can 
be seen that the projectile penetrated approximately 
1.27 cm (0.5 in.) . Figure 16 shows approximately 
2. 54 cm ( 1. 0 in. ) of penetration into a medium 
density (224 kg/m 3 ) target of aluminum foil and cork. 
Figure 17 shows a penetration of approximately 
5. p8 cm (2. 0 in. ) into a low density (72 kg/ m 3 ) target 
of aluminum foil and red polyurethane foam. The last 
target to be shown is in Figure 18. This is a very 
low density (32 kg/m 3 ) multilayer target. It is 
composed of 0. 25 mil doubly aluminized Mylar and 
red polyurethane foam spacers. It is the same 
material as that called the high performance multilayer 
insulation used in all the MDA configuration tests. 

This target was penetrated over 7. 62 cm (3. 0 in. ) 

It may seem intuitively obvious that a projectile 
would penetrate deeper into less dense target material. 
However, the important conclusions derived from 
these tests are In the calculations of the penetrated 
weight. That is, what weight target would be neces- 
sary to just stop a given projectile mass. This 
information is shown in Figure 19. The penetrated 
thickness is the black line with the scale on the left. 

The areal density or the weight per unit area 
necessary to stop the projectile is the dashed line 
with the scale on the right. Both parameters are 
plotted versus actual density of material. Looking at 
the areal density, it can be seen that the very low 
density (32 kg/m 3 ) high performance multilayer 
insulation would need half the weight that the medium 
and high density materials would need to stop the same 
projectile. This then is further evidence of the 
great effectiveness of high performance multilayer 
insulation in meteoroid protection. 

CURRENT RESEARCH 


In addition to doing penetration research on basic 
structures, the current research in the Materials 
Division Meteoroid Simulation Facility is concerned 
with studying the reactions of complicated structures 
to a hypervelocity impact. One method being used 
is to photograph the target during impact with a 
sequence camera. The camera being used is a six 
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Figure 14. Window after cleaning with alcohol. 
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Figure 15. Multilayers of aluminum foil and high 
density foam. 

frame image converter camera with a framing rate 
of up to 100 million frames per second and exposure 
times down to 5 nsec. 

Figure 20 shows one of the first series of 
pictures taken with this camera. The four frames 
were taken at 5-psec intervals. The picture at the 
upper left shows the projectile just coming into view 
in the grid. This grid was for measurement purposes 
and was positioned just to the front of a bumper. The 
second frame at the upper right shows the projectile 
just having penetrated the bumper. The expanding 
cloud of debris is already visible. The third frame 
at the lower left clearly shows the forward traveling 
debris plus the rearward traveling spray generated 
in the bumper impact. The fourth frame shows that 
the debris is not disappearing rapidly but must, in 
fact, be composed of some relatively slow debris. 


improved. Figure 21 shows two views of a bumper 
plus backup plate impact with 30 fisec separating 
the views. The picture on the left shows the target 
just after the projectile penetrated the bumper. The 
grid is behind the backup plate in this arrangement. 
The debris cloud is alroady expanding. The bright 
light on the outside is believed to be the long duration 
part of the impact flash. The view on the right shows 
the reaction of the backup plate. The filamentary 
debris is spall from the back of the backup plate. The 
tight nucleus of debris in the center is the small plug 
of material that was punched out of the backup plate. 

The optical camera is just one instrument with 
which a target impact can be studied. Other instru- 
ments that will be used are pulsed x-rays, strain 
gages, and accelerometers. The x-rays will be used 
for studying the inside of debris clouds and insulation. 
The strain gages and accelerometers will be used to 
measure the reaction rates of different structures. 
This is important because the strength of many 
materials are strain rate dependent. It is hoped that 
this type of information will make it much easier for 
structural design engineers to design meteoroid 
protection for spacecraft. 



'*j 


Since the foregoing pictures were taken, the 

optical arrangement for the camera has been greatly Figure 16. Multilayers of aluminum foil and cork. 
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Figure 17. Multilayers of aluminum foil and red 
polyurethane foam. 


CONCLUSIONS 

The penetration threshold for the basic MDA 
configuration is 0. 190 grams at 5.9 km/ sec. The 
threshold curve for this configuration is much less 
projectile -velocity dependent than is the curve for 
an equivalent single sheet. The introduction of a 
layer of high performance multilayer insulation into 
the MDA configuration approximately doubles its 
protective ability relative to projectile mass. 

The new sabot stripper design has been a 
great benefit to the penetration threshold programs 
by reducing ambiguous test results to a negligible 
amount. The improvements in launch techniques 
that allow more accurate velocity predictions have 
reduced the number of tests needed to arrive at valid 
penetration threshold data. 



Figure 18. Multilayers of aluminized Mylar and red 
polyurethane foam. 


Photographic and strain rate monitoring instru- 
mentation can be used to provide information on 
target reaction to impact. With this information, 
more effective meteoroid protection structures can 
be designed. 

Tests on various densities of multilayer 
insulations show that the very low density 0. 25 
mil doubly aluminized Mylar with red polyurethane 
foam is approximately twice as efficient as other 
medium density multilayer structures. 

The mechanical damage done to electrical cable 
protected by a single bumper is not excessive. 
However, the transient effects such as electrical 
arcing because of the plasma generated during 
impact, suggest that a better protection should be 
found. 
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Figure 21. Sequential views of a complete penetration of bumper and backup plate. 
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MHD ACCELERATOR RESEARCH 
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Patrick N. Espy 


A continual effort has been made in hypervelocity 
research to obtain ever-increasing particle velocities 
for impact testing. Because of their high potential 
flow rate, hydrogen and helium gases have been used 
to drive test materials. Using these gases, velocities 
of some 12 km/sec have been achieved for light gas 
guns. 

As higher velocities were sought, plasmas were 
considered. Hydrogen, air, aluminum, and other 
plasmas have been accelerated to high velocities for 
various types of testing. Included in these tests were 
hypersonic flow, thermonuclear research, and 
hypervelocity impact studies. Such plasma velocities 
were obtained through conical theta- and Z -pinch 
guns, rail, and coaxial J XB guns. 

Marshall Space Flight Center first entered the 
area of plasma-type hypervelocity particle accel- 
erators through a contract to MB Associates of San 
Ramon, California. As a result of the initial and 
subsequent fundings a particle accelerator was 
developed that produced particle velocities of as much 
as 15 km/sec from plasmas of 30 km/sec. A system 
of this type was delivered and placed in operation at 
Marshall Space Flight Center during this year. 

Figures 1 and 2 show the range and instrumen- 
tation rack of the original system. Figure 3 shows 
the capacitor bank. This bank consists of four 
modules that may be triggered simultaneously through 
four spark gap switches. A total of 240 kJ is 
available. These capacitors are dumped into the 


gun, which consists of two parallel copper rails 
shunted by a piece of aluminum foil. The foil is 
instantly vaporized, and the resulting plasma is 
accelerated by Lorentz forces of the current discharge 
through the circuit discharge and through the circuit 
formed by the rails and plasma. Small glass beads 
suspended between the rails are drag -accelerated by 
the plasma to a substantial fraction of the plasma 
veloc ity . 

After several months of operation in this mode, 
several. improvements became evident. The parallel 
rails were replaced by a pair of coaxial cylinders as 
shown in Figure 4. Only a small fraction of the 
capacitor bank was used initially, with careful 
attention given to matching the bank to the gun 
configuration. Velocities were obtained with this 
configuration that were comparable to those produced 
by the original device using practically full power. 

By a careful matching procedure, the bank was 
gradually increased until plasma velocities exceeding 
200 km/sec were measured. Even at these velocities, 
less than one-fourth the total available energy is 
being used. Figure 5 shows the range in its present 
configuration. Figure 6 is a photograph of the plasma 
emerging from the muzzle taken at an exposure of 
5 nsec with the B&W image converter camera. 

Having demonstrated our capability for producing 
extremely high plasma velocities, our research 
effort is now turning to better plasma diagnostics to 
determine plasma densities and temperatures, and 
to using this plasma to accelerate particles. 
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Figure i. Plasma range showing firing console. 



Figure 2. Plasma range from other side. 
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Figure 6. Self-luminous plasma emerging from coaxial gun. 
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IMPACT CODE DEVELOPMENT 

By 

R. J. Naumann 


Despite the recent improvements in accelerator 
performance, the prospects for development of a 
device that can launch fragile, low density meteoroid- 
like material at 72 km/sec are dim. Therefore, 
extrapolations of experimental data will probably be 
relied upon always. Such extrapolations are far 
more meaningful if they are based on theoretical 
results. 

Various computer techniques for treating the 
problem of hyper velocity impact are available. These 
have provided the theoretical basis for extrapolating 
to high velocities and have given a valuable insight 
into the processes involved in an impact. Bjork 
pioneered this effort in 1959 utilizing techniques 
developed by Landshoff and Harlow at Los Alamos. 
Since that time, various improvements have been 
made in equations of state and in putting more physics 
into the code. 

Basically, the code solves the inviscid hydro- 
dynamic equations in two dimensions. This requires 
the simultaneous solution of five nonlinear partial 
differential equations plus an equation of state in the 
presence of discontinuous boundary conditions. A 
set of mass plots for a particular case is shown in 
Figure 1. 

Bjork' s earliest work had a disadvantage in the 
fact that the determination of where the crater would 
stop growing was based upon subjective judgement. 
Later, under an MSFC contract, he improved this 
code to accommodate dissimilar materials and 
introduced a melting criterion to define crater depth. 

A series of runs was made of different material 
combinations at 20 km/sec and at 72 km/sec. 

In an attempt to avoid the large computer costs 
involved in the contract, Brown Engineering Com- 
pany as a support contractor was given the task of 
setting up this code for MSFC computers. A case 
was run for porous aluminum (p = 0.44 gm/ cm 3 ) 
impacting aluminum at 40 km/sec. The dynamic 
pressure profiles at various times are shown in 
Figure 2. Also shown is the residual material 
strength, which becomes zero at the point where 
sufficient irreversible work has been done to cause 


melting. In accordance with Bjork, this was taken 
as the crater depth . The results are shown in 
Figure 3. The point computed by Young at Brown 


(A) MASS POSITIONS ( B ) MASS POSITIONS ( C ) MASS POSITIONS 
t = 1 2m »*c + * 1 .9 u *«e + = A I m »•£ 
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♦ 7 A IS flii »*e * • -15 I i> 



CASE 8056. ALUMINUM INTO IRON AT 77 KM SEC 


Figure 1. Computer studies of impact phenomena 
represent the material by a series of imaginary 
mass points. As the impact progresses, the mo- 
tion of these mass points indicates the material 
flow. The regions where the mass points clus- 
ter to form the darker areas are regions in which 
the material has been compressed by the strong 
shock waves generated by the impact. Such shock 
pressures can be as high as 100 million atmos- 
pheres (many times the pressure at the center 
of the earth) . Material behavior at such pres- 
sures can be determined from quantum me- 
chanical calculations and has been measured at 
more modest pressures (1 to 2 million atmos- 
pheres) by controlled shock experiments. 
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Figure 2. Dynamic pressure profiles. 


Engineering Company agrees very favorably with 
extrapolations between Bjork's 11 km/sec and 
72 km/sec result. 

These types of codes do not take the strength of 
material into account in their solution, but as was 
seen in the preceding argument, the strength was 
entered in an after-the-fact fashion. Since the 
melting point was used to define the crater depth, 
such codes are automatically restricted to cases 
where strength is ignorable and where the crater 
depth is determined by melting. These restrictions 
precluded the check of such codes at available 
laboratory velocities. A strength code was developed 
by Bjork that predicted the crater depth of an 


aluminum — aluminum impact at 7.35 km/sec. It 
closely agreed with such a crater produced in the 
laboratory. Another interesting result from this 
effort was the prediction of the shock pressure as 
a function of depth in the target. Figure 4 compares 
Bjork f s results with measured pressures obtained 
earlier by GMDRL. Note the discrepancy at lower 
pressures. This is because of the method used at 
GMDRL to measure the pressure. The technique 
used was the throw-off-pellet technique developed by 
Rinehart. In this technique, the pellet traps the 
momentum imparted as the compression wave engulfs 
it. By measuring the velocity with which the pellet 
flies off, the shock pressure can be computed. 

Figure 5 shows the reason for the discrepancy. 
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Figure 3. Penetration versus impact velocity. 



Figure 4. Decay of the maximum value of <r at 
the axis with distance into the target. 


Characteristic of strength effects is the fact that the 
structure at lower pressures consists of an elastic 
precursor wave followed by the plastic wave. Appar- 
ently the pellet was trapping the elastic wave, thus 
flying off with a velocity that was too low. 

Another disadvantage of code work is the fact 
that the code results scale linearly with size, while 
actual craters do not. A slight tendency for smaller 
projectiles to crater less effectively has been 
observed in recent years. Considerable effort has 
been expended to determine the reason for this 
departure from size scaling. It was recently shown 


by Rosenblatt at Shock Hydro that strain rate effects 
could account for this phenomenon. Unfortunately, 
strain-rated data are not available for the strain 
rates involved in hypervelocity impact. When 
interpolating existing strain-rate data linearly, 
Rosenblatt obtained so much correction that the 
projectile actually bounced off the target. A 
logarithmic interpolation provided approximately the 
right correction, and Figure 6 shows a scaled large 
and small crater at two different times. The craters 
agree favorably with actual impact craters produced 
by particles of the sizes considered. 
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Figure 5. Two-wave shock structure approaching pellet at rear surface of target. 



z/D z/D 



Figure 6. Comparison of crater and splash for different size profiles 

at scaled times of U t/D =3.6 and 10. 8. 
o 
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SHOCK PROPAGATION IN SOLIDS 
IN HYPERVEIOCITY IMPACT 


By 

Tsun-Sen Fa 


INTRODUCTION 

To protect a space vehicle against meteoroid 
hazards, it is necessary to understand the meteoroid 
impact phenomena and the penetration mechanism. 
The meteoric particles, in spite of usually being 
small in size, possess high kinetic energies because 
of their hypervelocity nature (ranging from 
11 km/sec to 72 km/sec) . The term "hypervelocity 
impact" is defined generally as the collision of two 
solid bodies with a relative velocity higher than the 
elastic wave propagation speed in the material under 
consideration. 

The impact mechanisms are governed by an 
extremely complex system of nonlinear, partial 
differential equations with very complicated boundary 
conditions. Solving this problem analytically is a 
rather formidable task. Since 1958, several 
computer programs [1-6] have been developed that 
permit complete solutions through direct numerical 
analysis. A survey of this subject reveals that not 
only is a tremendous amount of computing time 
required for operating these computer programs, 
but large errors may develop because of the 
numerical fluctuations inherent in such solutions. 
However, the numerical approach is still considered 
to be the principal approach to obtaining a complete 
solution of hypervelocity impact. 

In addition to the numerical solutions obtained 
by the computer programs, analytic solutions of 
hypervelocity impact have been obtained by many 
authors. An excellent review of such solutions has 
been given by Rae [7] . So-called analytic treatments 
of hypervelocity impact usually refer to those 
approximate solutions that, even though less exact 
than the numerical methods, are easily obtained and, 
most important, still display the essential features 
of the problem. Evidently each approximate method 
is valid only in a certain regime, since solutions 
cannot be obtained in a rigorous manner by the avail- 
able analytic techniques without considerable 
simplifications. 


Blast wave theory [7,8] for instance, gives 
results in good agreement with direct numerical 
solutions at a late stage of impact. The blast wave 
solutions, however, are not applicable to the thin 
target plate since the thickness of the target is 
comparable to the projectile dimensions, and the 
point source assumption is no longer valid. The 
blast wave solution, however, does reveal the 
important fact that there exists a close analogy 
between the shock propagations from explosions and 
those from hypervelocity impacts. 

In the present study, a cylindrical projectile 
impacting on a target of infinite depth is considered. 
Immediately after impact, a plane shock is generated 
on the target surface; the shock area is equal to the 
circular cross section of the incident cylinder. 
Attention is then centered on the shock profiles as 
time elapses. This initial shock condition is 
considered to be equivalent to a surface explosion of 
finite area on the target surface. In fact, the 
surface shock equivalence has also been used by 
Zeldovich [9], among other Soviet scientists, as 
the idealized model that describes the initial stage 
of impacts. He did not, however, actually calculate 
the shock trajectories as functions of time. 

The formulations of the shock profiles in the 
present study are based on a shock dynamic theory 
[10 and 11] in fluid mechanics. The two-dimensional 
method of characteristics is employed to calculate 
the successive positions of shock waves as a function 
of time. The equations of state for metals are used 
in the formulation, and the shock interaction with the 
free boundary is also considered. 

The shortcoming of this theoretical model is 
that no information on the fluid is given behind the 
shock wave, except the flow properties deduced from 
Hugoniot relations. Because that shock wave develops 
in a semispherical shape in a short time after 
impact, quasi -similar solutions of blast wave theory 
[ 12] are used to approximate the fluid functions in 
the region behind the shock wave. 
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,, ,-Tho object of this study is to calculate the shook 
prof lion ro suiting from hypo rvo loci ty impact and to 
use thoso rosuitB In .predicting tho impact phonomona 
such as crater formation and spallation fraoturos. 
The applications of shock propagations In solid 
media will be roportod later. 

A numerical example is given in this report 
and compared with the existing results of other 
investigations. The agreement in general is good. 

The author wishes to express his appreciation 
to Mr. Carl Young for his discussions of the problem 
and to Dr. G. R. Guinn for his guidance and 
encouragement during the course of this work. 


THEORETICAL MODEL OF IMPACT- GENERATED 
SHOCK PROPAGATION IN SOLID MEDIA 


The physical phenomena of a projectile impacting 
on a target are illustrated schematically in Figure 1, 
In which a cylindrical projectile traveling at a 
velocity u^ is assumed to strike a semi-infinite 

target. Immediately after the impact (Fig. la) there 
are two strong shock waves moving in opposite 
directions: one travels forward into the target and 
the other backward into the projectile. The pressure 
is so high In the material through which the shock has 
passed that the material strength is negligible 
compared with the shock pressure. Therefore the 
material behind the shock waves is considered to 
behave essentially as an inviscid, compressible 
fluid. As the shock wave advanced further into the 
target material, the plane portion of the shock wave 
Is gradually consumed by the rarefaction waves from 
the comers (Fig. lb) . The shock front is distorted 
and attenuated in strength and eventually approaches 
the shape of a hemisphere. 

For extended elapsed time after impact, the 
crater is formed in the target while the isolated shock 
wave continues to propagate into the target. Further 
advancement into the semi -infinite target will cause 
the shock wave to degenerate into the plastic waves 
and then into the elastic wave. 


\I 
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AL 




EJECTED 

^/PARTICLES 
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X/I77T//777 


SHOCK WAVE 


RAREFACTION 
WAVES 

(a) Brief Elapsed Time After Impact 



(b) Intermediate Elapsed Time After Impact 



Figure i. General features of the hypervelocity 
impact mechanism. 



Because of the shock geometries, the problem 
resembles a moving shock configuration of unsteady, 
supersonic flow. Illustrated in Figure 2 is a normal 
shock moving along a channel whose wall suddenly 
enlarged at a certain section. As a result of the 


Figure 2. Moving shock front in a suddenly 
enlarged channel. 

enlargement, the shock wave changes its shape as 
it propagates into an enlarged area with a diminishing 
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of the plane shock by the rarefaction wave issuing 
from the corner of the wall. The mathematical 
formulation of this problem will be given in the 
next two sections. 

Figure la also Illustrates the manner in which 
the material particles arc splashed upward from the 
free boundary of the target, and the interaction of 
the shock waves with the target surface. The inter- 
action between shock waves and the free surface 
should be considered as the boundary condition in 
the calculation of shock profiles. (See page 50. ) 


FORMULATION OF THE CHARACTERISTIC 
EQUATIONS 



SUCCESSIVE 
SHOCK 
POSITIONS 
t - CONSTANT 


RAY TUBE 


The mathematical formulation of the shock 
profile relationships is based on a shock dynamics 
theory suggested by Whitham [ 13] for an originally 
plane shock with a corner. The two-dimensional 
computational scheme is essentially the same as 
that of Skews [ 10] , while the formulation of the 
axisymmetric case follows a recent paper of Hayes 
[ 11] . It should be noted that the formulations of this 
section only apply to a region away from the free 
surface of the target. In the neighborhood of the 
free surface, shock waves and free surface inter- 
action are considered separately and are presented 
in a later section. 


Two-Dimensional Case 


In Figure 3, an orthogonal coordinate system 
(/3, t) is constructed to describe the motion of shock 
waves, where t corresponds to the successive 
instants of shock wave history and (3 is called a 
ray to describe a specific portion of the shock wave 
in the course of traveling. The basic assumption is 
that each small portion of the shock, which is 
bounded by a ray tube, may be treated independently 
as one-dimensional flow and obeys the shock 
propagation law. The shock propagation law is a 
relation between the ray tube area A and the shock 
velocity U 

A = A (U) , (1) 

which will be derived in a later section. 


Figure 3. Ray tube coordinates system. 



Figure 4. Geometric relation of equations (2) and 
(3) for the two-dimensional case. 

PQ = U dt 
PS = Ad/3 . 

Let 0(/3, t) be the angle between the ray and a 
fixed direction. Then the change of 0 from P to S 
is 


In the physical plane (Fig 4) , the arc lengths 

PQ and PS are defined as 


60 


QR - PS 
PQ 


± 

U 
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or 


do. _ 1 _ . — n 

up “ u ' at 


i Similarly, 


SU - PQ 1 
■* - “T* " A 


OU M 

»P 6t 


or , 


M + ± . 9U 

at A op 


'(2. a) 


(2.b) 


(3) 


Notice that the derivations above are obtained 
for purely geometrical reasons . A combination of 
equations (2) and (3) indicates the kinematic relation 
between U and A 

±( _L . 9A \ . JL (± JU \ = ft 

9t\ u at / op y a a /3 J 0 • (4 * 

When equation ( i) is combined with equation 
(4) , the shock positions are determined for all 
times. The kinematic relations may be written in 
terms of U as follows: 


do 


d 


A'(U) 

U 


i 


ay 

at 


- o 


- du 

at A(u) a/3 0 


(5) 

( 6 ) 


where G is the equivalent sound speed as follows 
-U 


C 


A(U) A'(U) 


<!>) 


Lot 1 be tiie Illomann invariants of 

characteristics C and C „ respectively. Then 
equation (8) gives 


C + ; 6 + w(U) = i on ^ 
+ dt 


= C 


C “ ; 0 + w(U) = i on 


M = 


dt 


and 


"<«> ■ / Tc 


( 10 ) 

( 11 ) 

( 12 ) 


Let m(a t p) be the angle between a characteristic 
and a ray (Fig. 5) . The characteristic lines have 
the direction 


dy 

dx 


" tan (0 im) 


(13) 


on the physical plane. The positive sign in 
equation ( 13) corresponds to left-running character- 
istics, and the negative sign corresponds to 


or equivalently, 


a / A j (u) au \ a / l au \ 

at \ u at / a/3 \A(u) * dp ) 


= o 

(7) 


where A’(U) = 

du 


The above equations are analogous to the equa- 
tions of nonlinear sound waves or two-dimensional 
Steady, supersonic flows; hence the existing theories 
of compressible fluid flows on these topics can be 
equivalently applied to the present problems. 

By the theory of characteristics, equations (5) 
and (6) yield 



Figure 5. Mach angle in physical plane. 
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right-running characteristics; m Is equivalent to 
tho Much angle in two-dimensional steady flows. 


Figure (I shows that tho Increment IK} along a 
Left- running characteristic can be related to dt and 
d/1 In the (fi f t) plune as follows: 


dx -- U dt cos (0 -h m)/ (cos m) 
dy = U dt sin (0 + m)/ (cos m) 

dx = A d/3 cos (0 + m)/ (sin m) 
dy = A d/3 sin (0+ m)/(sin m) 


and 


, A d£ 

tan m = 7 ' dt 


AC 

(J 


(14) 


(15) 


Similar relations can be found for right-running 
characteristics by simply changing 0 + m to 
0 - m in equation ( 14) . 

Figure 7 shows that a plane shock, starting at 
the edge of a solid medium enters a sudden enlarge- 
ment area at time t = 0. For t > 0, there are two 
central fans emitted from the corners, which cause 
the shock front to expand. The normal shock 
portion is therefore diminishing in the course of time. 
Before time t (defined in Fig. 7), the physical 

picture is a simple wave pattern. Therefore U, 0, 




Figure 7. Sketch of shock propagation in solid 

medium (t is the time when the first characteristic 
c 

line reaches the shock centerline, x-axis, ; the 

normal shock portion vanishes completely 

for t = t ) . 
c 

and m are invariants along a characteristic line. 
Direct Integration of equation ( 14) yields 

x/t = U cos (0 + m)/cos m + x^/t , 

y/t = U sin (0 + m)/cos m+ y /t , 

° ( 16 ) 

x/p = A cos ( 0 + m) / sin m + x ^/ p , 

and, 

y/p = a sin (0 + m)/cos m + y J ol , 

where (x , y ) is the initial point. Similarly, 
o o 

changing m to -m in equation ( 16) gives the 
integrated equation for right-running characteristics. 

When t > t , the simple wave theory ceases to 
c 

apply and the shock profiles can be predicted only by 
constructing the characteristic nets, in the same 
manner as a two-dimensional nozzle flow field. 


1 ►x 

Figure 6. Relation of (x, y) plane and (/3,t) plane. 


The simple wave solution of two-dimensional 
cases is only hypothetical, since the shock waves 
generated by the projectile of finite dimensions are 
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three-dimensional flows. The simplest case for ' 
three-dimensional flow is the case of axisymmetry. 

In the next subsection, similar treatment shall 
be given for axisymmetric impact problems. A 
finite difference numerical scheme is adopted : 
because of the nonexistence of simple wave solutions. 
To do so, a starting line should be chosen to construct 
the net of characteristic lines. 1 The solutions on the 
starting line are the simple wave solutions described 
in this subsection. 


Axisymmetric Cases 


Let 0(/3, t) bo the angle between the ray and the 
axis of symmetry. The change of 0 from P to S 
(Fig. 9) is 


66 

RQ - PS 

1_ 

a a 

PQ 

u 

at r 

de 

1 3A 

A 


a/3 

Ur at 

+ 7 • 

sin 6 = 


where the relation (U at)/ 8r = sin 0 is used and 
0 = 0 at the axis of symmetry. 


As mentioned previously, the characteristic 
equations for axisymmetric shock propagation can 
be formulated similar to two-dimensional cases, 
except that the basic coordinate system (/3, t) must 
be redefined. 

To construct the axisymmetric solutions, the 
basic coordinate system (/3, t) is defined such that 
the distance along a ray between the shock positions 
given by t and t + dt is U dt, and such that 
A d/3/ r is the distance between the ray (3 and 
/3 + d/3; r is the distance from the axis of symmetry 
and U is the shock velocity (Fig. 8) . 



Figure 8. Basic coordinate system of 
axisymmetric case. 



Figure 9. Geometric relation for 
axisymmetric cases. 


A similar treatment for the change of 0 from 
P to Q gives 


30 r 

at a 



(18) 


Let 61 = A6/3/r be the distance between the rays 
(3 and (3 + <5/3. Equations (17) and (18) th^n become 


de 1 dA sin 6 
dl ’ A * U 8t + r 


( 19 ) 


i. In steady, axisymmetric supersonic flow calculations at the exit of a nozzle, the Prandtl-Meyer 
expansion fan, which is a simple wave solution, is used to provide such a starting line f 14] . 
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and 


• Conservation of mass 


DO OU 

u at + u M 


o 


( 20 ) 


Assume that the shock propagation law, i.e. , 
the shock speed and area relation, is expressed in 
the form (see page 54). 


dA 

A 



^ m = M(U). 


Po u = p (U -u) (24) 

• Conservation of momentum 

P - Po = Po u u (25) 

• Conservation of energy 

e-e 0 = -i- (ft+P )(^ - 7 ). (26) 


Equation ( 19) , then, can be rewritten as 
90 1 9U sin 0 

da ~i? ' u 2 at r ° ■ 


where U, u, p, p, and e are the shock velocity, 
fluid velocity, density, pressure, and the specific 
internal energy of the fluid, respectively. The 
( 21 ) subscript 0 refers to the undisturbed (or stationary) 
region. 


Applying the theory of characteristics in 
equations ( 20 ) and ( 21 ) yields 

r±: ^ ±de + p — - U dt = 0 (22) 

pU r 


One more relation is needed to solve the flow 
across a shock; that is, the state equation of the 
material 

e = f (p, p) (27) 


along 


Cl : 


di 

U dt 




(23) 


where r± is the compatibility equation of the 
physical characteristics C-j=. 

For a given initial shock shape, the shock 
shape at later time can be obtained by constructing 
a finite difference net in the same manner as that 
of a steady supersonic jet downstream of the nozzle 
exit section. 


Recently, many studies have been conducted [15-17] 
on the state equation of solid material. So far, no 
simple and general equation has been formulated. 
For the sake of convenience, an equation, which 
relates the shock velocity and the particle velocity 
u, is used in this work. In general, there are two 
forms; one applies to so called c, s material and the 
other applies to aluminum. 


c, s Materials 


Two boundary conditions are imposed. One 
is simply that 

0=0 on r = 0 

The other condition is on the free boundary where 
the shock wave interacts the free boundary. The 
interaction will be described in a later section. 

HUG0NI0T RELATIONS OF SHOCK WAVES IN 
SOLID MEDIA 

When a shock wave moves into a stationary 
medium, the flow variables across a shock are 
governed by the following laws: 


In Reference 16, McQueen et al. use a linear 
equation, 

U = c + s x pu , (28) 

which relates the velocities U and u for most 
metals, c is dilatational wave speed and p is a 
constant that is tabulated for 16 materials in their 
paper. The materials that obey equation (28) are 
classified as c, s materials. 

Substituting equation (28) into equations (24) , 
(25) , and (26) , the following flow variables are 
obtained: 

U = (U - c)/s , (29) 
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S U 

p (s-l)U + c ^ 


(30) 


and 


U (U - c) 

P - — ^ Po 


(31) 


The slope of the Rankine-Hugoniot curve in the 
(P, p ) plane is given by 


C\ (2 M s - 1) [(s - 1) M + l] 2 , 


2 c 
c-H V 


(32) 


where M = U/c is the shock Mach number, 
s 

differs slightly from the adiabatic sound speed 

as M nears i, but appreciable discrepancy is 
s 

found for large shock Mach numbers. 

The adiabatic sound speed, a = (9p/c)p)s, 
should be derived by an isentropic process that 


requires the explicit form of the state equation, 
equation (27) . But for c, s material, the adiabatic 
sound speed, a, can be obtained approximately in 
terms of the material constants c and s. The 
resultant expression is 

e 2 

a 2 = - [(3s-i) (s-1) M g 2 - 2 (s 2 - 3s+ 1) M g 

- (2s - 1)] . (33) 

A detailed derivation is given in the appendix, 
and a comparison of a and c^ is made in 

Figure 10. 

It is of interest to point out that the c, s material 
is very similar to an ideal gas. The range of s 
varies from 1.2 to 1.9 for most metals, which is 
analogous to the specific heat ratio of gases. It 
should be stressed, however, that the c, s relation is 
only an empirical equation. Its accuracy for the 
shock pressure, greater than 2 mb, is highly 
questionable. 



0 2 4 6 8 10 12 14 16 18 20 22 

U/C 

Figure 10. The ratio of adiabatic sound speed and the slope of Rankine-Hugoniot curve as a function of 

shock Mach number of c, s material. 
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Aluminum 


Aluminum is one of the exceptional metals that 
do not follow the linear shock relation of equation 
(28) . Instead, the shock relation for aluminum Is 
expressed in a quadratic form as 

U - aj + a 2 u + a 3 u 2 (34) 

In addition, the adiabatic sound speed a Is also 
a quadratic function of u; that is, 

a = b t + b 2 u + b 3 u 2 . (35) 


Consider a shock wave moving in a channel with 
varying cross-sectional area. The shock strength 
will attenuate if the channel area increases down- 
stream; and conversely, the shock wave will bo 
strengthened if the channel area converges. This 
I Genome non is well-known in blast wave theory. 

Even if the fluid is assumed frictionless, the shock 
front velocity decreases very quickly and approaches 
the velocity of sound at a large distance from the 
source of explosion. The enlargement of the shock 
front area results in attenuation of the shock strength, 
or in the equivalent, a reduction of shock velocity. 

The shock propagation law relates the change of the 
shobk velocity to the change of the channel area. 


For simplicity, the coefficients a. and b. used 

in this study are obtained by fitting them into the 
table of Chou et al. [ 18] (Table 1) . 

TABLE 1. SHOCK AND ADIABATIC SOUND SPEED 


CONSTANTS FOR ALUMINUM 



u ^5.0 ( mm/psec ) 

u > 5.0 ( mm/psec) 

a i 

5.985 

5.278 

*2 

1.237 

1.503 


-0.001 

-0.0271 

b i 

6.883 

5. 101 

b 2 

0.9849 

1.518 

b 3 

-0.0084 

-0.048 


Chou’s table was established by using the 
Tillotsen equation [14] (for p > 1 mb) and the 
experimental data of the Los Alamos Laboratory 
[ 16] ( for p ^ 1 mb) . 

Again, equation (34) , combining equations (24) , 
(25) , and (26) , provides the shock relations needed 
in the next section. 


The shock propagation law was first obtained 
by Chester [19] and Chisnell [20] from a linearized 
approximation of one -dimensional compressible 
flows. Their results were later rederived by 
Whitham [21] who applied a simple rule in the non- 
steady, one-dimensional flow theory. Their result 
is sometimes called CCW Approximation. In this 
section a shock propagation law is formulated for 
solid media. 

Whitman’s rule can be stated as follows: the 
flow quantities along a characteristic line must be 
equal to the flow quantities just behind the shock 
wave. With this rule and the Hugoniot relation, a 
differential relation is obtained that describes the 
motion of a shock wave down a nonuniform area 
channel. In hypervelocity impact problems, a 
plane shock is generated initially at the interface 
between the projectile and the target, and the shock 
front area suddenly enlarges while it propagates into 
the target. Then the shock propagation law is applied 
under the assumption that the whole flow field is 
composed of a number of one-dimensional rays 
(see page 43), each one independently obeying this 
law. 


For one -dimensional, unsteady flow in a 
variable area channel, the governing equations of 
the fluid flow may be written in the following form: 


at ax p \ax a dx / 


o 


(36) 


SHOCK PROPAGATION LAWS 


As mentioned previously, a shock velocity and 
area relation (equation 1) or shock propagation law, 
is needed for completing the characteristic 
solutions. Equation ( 1) is derived in this section. 


da , Bu a 2 Dp 

— + u — + — * — tL 

at dx p dx 


0 


(37) 


and 


Be 

at 


+ 




+ u 


dx 


) 


0 . (38) 
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Applying the theory of characteristics to 
equations (36) and (37) yields 




du + A 5i + a) dA = 0 (39> 


For aluminum, the function p 2 of equation (41) is 
, U (u+a )l __1 


u 


U “ U-u (U-u) (a 2 +2a 3 u) 


on 


C 


dx 

dt 


= u ± a 


where denotes hodographic characteristics and 
C ± denotes physical characteristics. 

To apply Whitham's rule, the flow quantities on 
the outgoing characteristics C + are matched with 

the flow quantities behind the shock wave. The 
appropriate branch in equation (39) is 


+ 


1 

a (a 2 + 2a 3 u) 


(44) 


where a and u are functions of U, as shown in 
equations (34) and (35) of the previous section. 

Integration of equation (40) yields, K being 
an integration, 

A(U) = Kexp- / M 2 ^ (4 

which is equation ( 1) on page 43. 


df> + a.dfl + du + du + d£ = 0 

p u p U a A 

The terms p, a, and u are functions of 
shock velocity U and the material characteristics. 
Equation (40) can thus be integrated, and the 
integrated equation represents the shock propagation 
law. 

For c, a material, equation (40) becomes 


Notice that the integration in equation (45) can be 
carried out explicitly through a tedious procedure. 

In actual computations, p is used instead of A(U) . 
The computational procedures will be discussed in a 
later section. 


SHOCK WAVE INTERACTION WITH 
FREE BOUNDARY 


T " - * (M s> IT ’ (41) 

where M = U/c and 
s 

1 = 1 + s 

7 ? [ i+ (s - 1) M ] M -1 

s s 

M 

♦ 2 X 

{[(s-l)M +1] [(3s- 1) M - (2s -l)]} 2 

S o 

{[(s-i)M g +i] [ (3s - i) M g - (2s - 1) ] } i 

+ (M -i) [ i+ (s - 1) M ] 

S S (42) 


As stated previously, information about shock- 
free surface interaction is needed to specify the 
boundary condition in constructing the characteristic 
solutions. Figure 11a shows an expansion fan behind 
the shock, through which the pressure drops from 
the value immediately behind the shock wave to zero 
at the tail of the fan. The physical picture is analogous 
to Prandtl-Meyer expansion around a convex corner 
except that the expansion fan is three-dimensional and 
unsteady in nature. The shock front would attenuate 
and is distorted by the rarefaction waves. 

A simplified model has been proposed by 
Russian scientists [22] to estimate the shock 
inclination. Recently, Rae [7] applied it to the 
impact-generated shock propagation, but he chose 
ideal gases in the formulation. In the present study, 
Rae's formulation is generalized to solid media. 


For the case of strong shock, M = U s/c — 00 , 
equation (42) reduces to 




2 

00 


[ ( s - 1) (3s- 1)] * 
1+ [(s-i) (3s- 1)]* 


(43) 


From Figure lib, if the shock inclination (3 is 
greater than a certain value 0* called the critical 
angle, the rarefaction waves cannot reach the shock 
front. Therefore, if the original shock inclination 
13 > (3*, the shock motion along the free surface will 
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SHOCK l ROUT 


(a) EXPANSION BEHIND THE SHOCK FRONT 


and 


-1 

/a 2 

2u 

, U 2 \ - ] 

COS 

(li 7 

U 

' uV J 


Equation (4G) indicates the critical shock 
inclination angle ft* 9 as a function of U, u, and a. 
Since a and u can be expressed in terms of 
U, ft* actually depends on shock velocity U only. 



AT t 

(b) CRITICAL ANGLE B* 

Figure 11. Shock wave and free surface interaction. 

not be affected by the free surface interaction. 
Conversely, if ft > ft*, ft will change continuously 
and eventually approaches (3*. The shock will 
remain undistorted afterward. 

The transition phenomenon (ft— ft*) is rather 
complex in mathematical analysis, because an 
unsteady and three-dimensional flow problem is 
involved. In this work, it is assumed that the 
transition period is so short in the development of 
shock profiles that it can be neglected. Then 
ft - (3 * will provide the boundary condition at the 
free surface, ft* depends on the local shock strength 
and is derived below. 

Referring to Figure lib, the cosine law gives 

BC 2 = AB 2 + AC 2 - 2 AB • BC cos 0* 
or 

a 2 = u 2 + U 2 / cos 2 p* - 2uU 


CONSTRUCTION OF NUMERICAL SOLUTIONS 
FOR SHOCKWAVE PROPAGATIONS 


Initial Normal Shock Waves 

When a projectile impacts on a target with 
relative velocity V, two normal shook waves arc 
generated at the interface and move in opposite 
directions with respect to the interface (Fig. 12a). 
The normal shock solutions can be obtained easily by 
applying equations (24) , (26) , and (28) or (34) . 

For similar material hypervelocity impacts, the 
particle velocity equals one -half the impact velocity 
(Fig. 12b), 

u 3 = u 2 = V/2 . (47) 

Equation (47) is derived by taking a vciocity 
transformation. If the projectile and the target move 
in opposite directions with the same speed V/2, 
the fluid particles behind the shock waves should be 
motionless, because of symmetry, with respect to 
the laboratory coordinate system. 

For nonsimilar material impacts, a contact 
discontinuity interface exists. The pressure and 
velocity functions are assumed continuous, however, 
across the interface; P3 = P2 and u 3 = u 2 . From 
equation (2b) 

P2 = Po u 2 u o (48) 

and 

Ps = Pi (V-u 3 ) (V-U t ) 
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(a) NORMAL SHOCKS GENERATED IMMEDIATELY AFTER IMPACT 


PROJECTILE 

" 1 " 


uj = V 

Pj = 0 
Pi 



TARGET 

"O" 


" 3 " 





Figure 13. Construction of a simple wave solution. 


0 X is 



dU 

Uji (U) 


(30) 


The shock position of (joint 1 is determined by 


(b) NORMAL SHOCK MECHANISM 


t 


Ui cos (0i + mjl/cos mj 


Figure 12. Impact-generated normal shockwaves. 

With the aid of shock relations, equations (28) 
or (34) , the solutions of V x and U 0 are readily 
found . 

Simple Wave Solutions 

As shown in Figure 13, U 0 is assumed to be the 
normal shock velocity claculated in equation (47) or 
equation (48) . Then the Mach angle m 0 can be 
obtained from equations ( 15) and (41) , 

mo = tan 1 [p (U 0 ) ] , (49) 

and 


and 

sin (0j + m^/cos mj , (51) 

where t is the time after impact. 

Repeating the procedure, by choosing a 
sequence of shock velocities U 2 , U 3 , U^, m., tf., 

x /t, y /t (i = 2, * * * * n) are determined in the 
i i 

same manner. The computations continue until 
x./t = 0. This completes the construction of a 

simple-wave solution. 

Asymptotic Solution of Two-Dimensional 
Shock Profiles For Strong Shock 


o 0 = 0 

An analytic expression of shock profiles can be 
Next, choose a shock velocity Uj, Uj < U„, obtained for large U in simple wave solutions, c, 

and calculate m 0 in equation (49) . The shock angle s material is chosen in the formulation that follows. 
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The function p(U) for a c, s material in the 
shock propagation law approaches a constant as 

(equation 43) a and is denoted by Thoro- 

Ibro, 


..j- r f ■ 


•7- *n (U/U,) 

'on 


and 


H 0 

' on 


U = U 0 e 


(52) 


The shock wave locus is in the form 
/! 6 

x 00 / 

—7 = e cos (0+m ) /cos m 
U 0 1 00 00 


and 


U 0 t 


fJL 6 

00 / 
e sin ( 0 + m )/cos m 

oo' 00 


(53) 


where m = tan 1 a is also a constant. 

OO ^OO 

Lotting R= (x 2 +y 2 ) 2 , equation (53) maybe 
expressed in polar coordinates as 


R - 




k e 


(54) 



SHOCK 

PROFILE 


x 


Figure 14 . Approximate shock profile for the strong 
shock situation ( U large) . 


• 0 - 0 at the axis of symmetry 

• 0 = - ^ p*+ on the free surface. 

The method of characteristics is applied in 
computing the shock profiles. To do this, the 
characteristic equations, equations (22) and (23), 
are first written in finite difference forms, and the 
solutions can be constructed by the standard technique 
in the compressible flow theory [ 14] . 


where 

(p = 0 + m 

^ oo 

and 


cos m 

oo 

Equation (54) (Fig. 14) represents a 
logarithmic spiral curve, and has been derived by 
Skew [ 10] for an ideal gas. However, it should be 
emphasized that equation (52) , must be used with 
caution since the simple linear shock relation of 
c, s material may not be valid for very strong shocks. 


The finite difference scheme developed herein 
differs from the conventional method. The present 
scheme is designed to construct the shock profiles 
for successive time intervals. Instead of using two 
known field points to obtain the third new point, only 
one point and the time interval At are required to 
determine a new point. The advantage in this 
numerical scheme is that the successive shock 
profiles as a function of time can be directly com- 
puted without interpolations between the field points. 

The computations for this study are programmed 
on the IBM 1130 digital computer. 


Finite Difference Solutions of the 
Characteristic Equations 


For an axi symmetric impact problem, the 
simple wave solution only provides the initial con- 
dition. Two boundary conditions are used: 


AN APPROXIMATION FOR PRESSURE PULSE 
AT LATE STAGE OF IMPACT 


The shortcoming on the present shock 
propagation model is that no solutions about the 
flow variables are given behind the shock wave. The 
flow field of a diffracting shock wave is too 
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complicated to be obtained analytically. This Ik tlio 
cHHontlal roaHon for the development of direct 
numerical progrumH, Hindi an PIC code and 
Kuloroan c(xlo. In fact, the flow variables aro 
significant in the description of the impact mechanism. 

At the late stage of shock propagation, it is 
found that the shock front develops rapidly into a 
sphere-like shape, and the flow variables vary 
slowly along the shock front. These facts suggest 
that it is possible to obtain an approximation for the 
flow field at the large elapsed time. 



To find the approximate solutions in the flow 
field, the equations for three-dimensional, unsteady 
flow with symmetric axes are written in a spherical 
coordinate system. These equations are 


Figure 15. Spherical coordinate system with 
symmetric axis. 

By substituting equation (56) into equation (55), 
the zero order equations are obtained as follow: 
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(57) 


+ — x cos (p - 0 , 
r 


(55) 


and 


In the above system of equations, u = u^ and the 
subscript o is dropped. 


de 

at 


+ u • — 


de 

ar 


u </> . ae 
r dtp 


+ . is 

r a <p 




M 

at 


+ u 


ifi 

r at 


Equation (57) is identified as a one-dimen- 
sional, unsteady flow equation system. A trans- 
formation of the dependent and independent variables 
is introduced as 


where u^ and u^ represent the velocity com- 
ponents in the r and <p directions, respectively 
(Fig. 15). 

Expanding the variable u^, u^, p, and P in 
Fourier series yields 


u r (r, <p, t) = u q (r, t) + Y, “ (r, t) cos n 0 , 

n=i (56) 


u 0 (r, cp, t) = £ 
n=l 


v (r, t) sin n 0, etc. 
n 


u = U £ 


P = P Q V* g . 
p = P Q h , 


« "a * 

, _ R dU 
U dR 


and 




(58) 
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whore II Is the shock front position. The first two 
equations of equation (57) become 


It t\ 1. hi£ i - X II — 

(f '*> HI +h H { ~~ xv w 


(59) 


(f 


‘>5 


+ A£ *if * - 


-I 


(60) 


Note that f, g, and h are functions of two 
variables, 4 and R(t). To solve the above system 
of equations is still a very difficult task. In general, 
the similarity solutions do not exist because the 
assumption of a strong condition may not be justified. 
In fact, when the impact-generated shock wave 
becomes spherical in shape, the shock strength is 
usually weak. 


Now, two assumptions are imposed as follow: 


In terms of physical quantities, the solutions are 
u = uj r/lt , 
p - pi(r/ll) m 


P 


where 

m 



3u,/U+\U/ Pl - 
(ui - U)/U 


and 


a = * 

p* 


The subscript i denotes the Hugoniot flow quanitities. 


f = M 

ah = Jh_ _ dJhj 

au h t * du 

where f lf h 4 are the functions evaluated at the 
shock front (4^=1). They are time-dependent 
functions. 

The first assumption is made because the 
numerical solutions of other studies [8,23] show 
that the velocity functions are nearly linear 2 , while 
the second assumption is based on a quasi-similarity 
theory [12]. 

From equations (59) and (60), combined with 
the assumptions above, g and h are readily 
obtained as follow: 

g = A (| m+2 - l)+ gi 

and (61) 

h = h t | m 

where gj = g( 1, t) and A, m are functions of 
ft, gi, hj, and A. 


NUMERICAL RESULTS AND DISCUSSION 

In this section, some numerical results are 
presented. Based on these results, discussions are 
given on the theoretical model proposed in this study. 
Since some of the assumptions imposed in this model 
are subject to a certain degree of uncertainty, it is 
necessary to compare the present results carefully 
with numerical data available from other studies. 
Because of the complexity inherent in this problem, 
it is believed that the justification of some of the 
assumptions is possible only through comparison. 

The Starting Point of Shock Curvature 

In Whitham's rule, the effect of disturbances 
behind the shock wave is ignored. The justification of 
this model is difficult to establish, because the 
information on the flow field behind the shock wave is 
still meager. 

Point A in Figure 16a represents the location 
where the shock curvature starts. By the CCW 
approximation, the following is obtained: 

tan m 0 = p(U 0 ) (62) 


2. The first assumption was also made by Rayzer [24] and Sakura [12] , but no basic reasons were given. 
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But Mq can also be predicted by the following 
alternative equation: 


tan m 0 

:.'J i 



(63) 


Equation (63) was derived by Skew [25] and also 
by Heyda [26] to indicate the point at which the 
normal shock wave is first reached by the rare- 
faction wave from a turning corner. It must be 
pointed out that equation (63) is only an approxi- 
mation because the shock wave is not uniform and 
the velocity u and sound velocity a depend on the 
locations. 


i 



A comparison of the above two approximations 
is shown in Figure 16b. For high impact velocities, 
the discrepancy between these two predictions is 
small, but a large discrepancy exists in the low 
velocity region. Skew [25] performed an experi- 
ment for shock waves diffracting around a corner. 
Observation of the shock shape favors the prediction 
of equation (63) . This may imply that the CCW 
approximation only works for high impact velocities. 
However, this is not a definite conclusion since the 
phenomenon of impact-generated shock waves may 
be different from Skew's experiments. 


(a) 



The Successive Positions of 
Shock Profiles 

A numerical example of aluminum impact on 
aluminum at a velocity of 20 km/sec is presented 
in Figure 17 such that the positions of shock profiles 
are functions of successive time steps. The normal 
shock portion vanishes very rapidly, and it completely 
disappears at about 5. 5 psec. The shock profile 
develops further into an ellipsoid and eventually 
approaches a hemispherical shape. 

The appearance of hemispherical shock shape 
inspires the assumption of spherical symmetry 
made by many investigators (7) . This assumption 
is reasonable only if, after a long period of time, 
the shock strength has attenuated considerably and 
degenerated into the speed of sound. Soon after 
impact, the shock shape is far from spherical in 
shape and, most important, the flow properties 
vary along the shock front. The pressure distri- 
bution, for instance, has the highest value at the 
symmetric axis and decreases in value toward the 
target free surface. 


10 h 


0 1 — — - I 1 I [ . | 

0 20 40 60 80 100 

IMPACT VELOCITY (km/sec) 


Figure 16. The starting point of shock curvature 
in aluminum. 


Pressure Distributions 


A comparison is made with the numerical results 
obtained by Heyda and Riney [27] concerning the peak 
pressure along the axis of cylindrical symmetry. 
Figure 18 shows that for (z/d) < 1, the pressure is 
equal to Hugoniot pressure obtained by the normal 
shock relation since the plane shock generated at the 
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interface has not been attenuated by the rarefaction 
waves from the corner. The pressure then decreases 
drastically as the shock wave advances into the 
targe l . The discrepancy of the so two results at the 
early stage is mnjniy a result of tho numerical 
diffusion effect inherent in tho direct numerical 
scheme. Krrors of about 90 percent occur compared 
with the exact solutions. However, tho agreement In 
general is good except the late stage In which the 
present solution is consistently higher than the direct 
numerical solutions. This is believed to be caused 


by the fact that the rarefaction waves reflected from 
the rear surface of the projectile are not considered. 
Such rarefaction waves are known to attenuate the 
target shock strength [ 18] . 

Tho angular distributions of shock pressure 
are shown in Figure 19 so that the pressure gradients 
in tho <p direction are indeed small for largo 
elapsed time. This phenomenon provides the basis 
of formulating the quasi-similarity solutions in the 
previous section. 


R/d 



Figure 17. Numerical example of aluminum impact on aluminum at a velocity of 20 km/sec. 
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APPENDIX: SOUND SPEED OFc, s MATERIAL BEHIND THE SHOCK 


The Mie-Griineisen equation of state is 

P - Po = P r (p) (e - o c ) , (A-l) 

where the subscript o denotes the initial state and 
l'(p) , the Gruneiscn ratio, is a function only of 
density. Differentiating equation (A-l) with respect 
to c at constant p, the following expression for 
T(p) in terms of thermodynamic quantities is 
obtained: 


r 


t 

p 




(A-2) 


According to Reference 7, r « 2 s - i at 
normal density, or pressure p = 0. But for higher 
pressure, a constant value of T was proposed by 
Rao [7], who matched the Mie-Gruneisen equation 
of state with the Rankine-Hugoniot relation of c, 
s material at very high shock strength. It is 


T -- 2 ( s - 1) 


( A-3) 


A general expression of the equation of state 
may be written 

P - p(p, S) ; 

hence , 



where [ ^ j — a is the adiabatic sound speed, and 

V 9 p/r 



( A-7) 


Substituting equation (A-3) into equation (A-7) 
yields 



pT r 


( A-8) 


It should be pointed out here that the range of 
validity of this expression of T is uncertain since 
there is no way to justify it by means of the available 
experimental data of metals. In spite of this, 
equation (A-3) is used in this work for the sake of 
convenience. 


Now, a thcr modynamic relation can be derived 
to yield the adiabatic sound speed a as a function 
of dilatational wave speed c and the material 
parameter s. The derivation presented here is 
based on Reference 28. To do this, the Rankine- 
Hugoniot equation is differentiated first as follows: 


de 


_1 

2 



dp - (p + 



(A-4) 


Then equation ( A-G) can be written as 


dp = a 2 dp + p r T dS 


(A-9) 


Eliminating the term T dS from equations 
(A-5) and (A-9) gives a relation between the 
adiabatic sound speed, a 2 , and the slope of the 
Rankine-Hugoniot curve in (p, p) plane, c 2 ^ , 




Combined with the second law of thermody- 
namics, equation (A-4) becomes 



(A- 10) 

The shock relations of c, s materials can be 
expressed as follows: 

U = c + sU 


Pi _ s M 

p 1 + (s - i) M 
o 


(A- 11) 
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Pl ^ = M (M - 1) 


pc s 


and 


2 _ dg _ dP/c/M 

C Ii-H d p " dp/c/M 


= (2M - 1) [1+ (s - 1) M] 2 , 


Substituting equation (A-ll) inequation (A-10) 
yields 

P on 1 

a = — [ (3s - 1) (s - 1) M 2 - (2s 2 - 0s+ 2) M - (2s - 1) ] a 

(A- 12) 


and 


"R-H 


(3s - 1) M - ( 2s - 1) 
(2M - 1) • i+ (s - 1) M 


-i 1 

2 


. (A- 13) 


where 1V1 - U/c . 


The ratio of a/c ^ depends upon s and M only. 
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ONE-DIMENSIONAL MODEL FOR INVESTIGATION 
OF CHANGES IN CONDUCTIVITY 
OF DIELECTRICS UNDER SHOCK CONDITIONS 

By 

E. E. Klingman 


ABSTRACT 

A one-dimensional model is used to investigate 
changes in conductivity of dielectrics under extreme 
shock conditions. The dielectric is forced into 
metallic conduction states by varying the energy 
band structure as a function of Lattice parameter. 

The most striking feature of the treatment is its 
generality, deriving from the fact that materials in 
which repulsive forces predominate can be described 
in much the same way. This is in contrast to the 
various bonds that must be considered when 
attractive forces arc most important. Although the 
model is qualitatively successful, it is concluded 
that a nonlinear treatment is necessary to obtain 
quantitative results. 

INTRODUCTION 


The Pegasus Meteoroid Detection Satellites use 
large surface area capacitors as detection devices. 

A meteoroid-hit upon one of the capacitors results in 
a current pulse that furnishes information to the 
telemetry system. The mechanism producing the 
pulse may be one of the following: 

• Penetration and physical shorting of 
capacitors 

• Plasma generation 

• Dielectric shocked into conduction 

When this study was begun, the possibility of the last 
mechanism was uncertain. The object of the study 
was to investigate the possibility of dielectrics 
being shocked into conduction and, if possible, to 
relate the pulse duration to meteor velocity. 


Experimental data [ 1, 2] were found that 
established beyond doubt the possibility of forcing 
insulators into conduction states. The application 
of sufficiently high pressures to a variety of 
insulators resulted in conductivities that approxi- 
mated those of metals. The study then centered 
on finding a suitable model for calculations. This 
model would then be used to investigate velocity 
effects. 


MODEL FOR CALCULATIONS 


The situation, even in the best of crystal 
structures, is formidable from the viewpoint of 
realistic calculations. It was quickly decided that 
any attempt to analyze the actual dielectric, Mylar, 
would soon bog down in complexity. In view of the 
time available for this study, a one-dimensional 
model was chosen with the hope of obtaining 
qualitative information that could be correlated with 
data from the satellites, and experiments already 
performed on the ground. Although symmetry plays 
a vital role in crystal calculations, the added 
complexity was of questionable value because of the 
unknown effects on the symmetry of the high 
temperature and pressures involved during collision. 

The differences in conductivity of metals and 
insulators vary over an extremely large range (on 
the order of 10 42 ) . This variation can only be 
accounted for in terms of the structures of energy 
bands. The simplest model resulting in band 
structure is the Kronig-Penny model [3], but this 
model was found inadequate for pressure calculations. 
Since pressures of 1 Mb compress condensed systems 
by a factor of 2 [ 1] , the parameter to be varied with 
increasing pressure is the lattice constant a. The 
increase in temperature is taken into account through 
an approximation that will be explained later. 
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Conditions are such in the case of extreme 
shock that considerations normally of importance 
must he evaluated. Of critical Importance in crys- 
tal calculations at normal pressures and temperatures 
Is the ty|)o of crystal involved; i.o., Ionic, molecular, 
covalent, metallic. The typo of crystal is deter- 
mined by the nature of the attractive forces causing 
the atom to crystallize. In shocked materials, 
however, energy exceeds the attractive energy, 
and the repulsive potential between the atoms predom- 
inantly determines the properties of the system. 
Empirical results from shock data indicate that at 
milder conditions than expected, it is permissible to 
neglect the details of binding of a substance f 4] . 

This could have been predicted from the success of 
the Fermi -Dirac model which ignores band structure 
and statistically groups atoms about the nucleus. 

When shock and heat energies approach the chemical 
binding energies, this model becomes applicable since 
at that point all materials become close-packed and 
metallic. 

At this point a digression may be in order con- 
cerning the difference in metals and insulators from 
an energy band standpoint. 

Bands arise from the fact that no two particles 
obeying Fermi statistics can be in the same quantum 
state in a given system. Consider, for example, 
three atoms in their ground or lowest energy state. 
Quantum mechanics predicts for each isolated atom 
a set of energy states that are identical to those of 
any similar isolated atom. The situation changes 
when these isolated atoms are brought together to 
form a system. As the wave functions overlap and 
interact, the states are shifted to avoid violating the 
Pauli exclusion principle. The group of states that 
were identical for similar isolated atoms now form 
a band. This is shown in Figure 1 for six atoms. 

In a crystal with very large numbers ( 10 23 ) of 
atoms, the levels are so closely spaced inside a band 
that the band can be considered continuous for all 
practical purposes. 

A crystal composed of N atoms is presented in 
Figure 2. The particular atom chosen has two 
electrons in the first energy state and one electron 
in the second state. The lowest energy band will 
have N states, each of which will be occupied by 
two electrons with opposite spins. The opposite spins 
prevent electrons with the same energy from having 
identical quantum numbers, thereby violating the 
Pauli exclusion principle. The 2N electrons fill the 
first band to capacity. The next band, with only N 



Figure 1. Formation of bands. 


CRYSTAL 

SPACING 



interatomic distance 


Figure 2, Energy bands for crystal with N atoms 
of is 2 2s* configuration. 

electrons and 2N available states, will be only 
half-filled. Electrons in this band can gain thermal 
energy and make transitions to higher states in the 
same band. If an electric field is applied to the 
crystal, conduction is possible in the half-filled band 
No conduction is possible in the filled band. 

If the atom chosen had been one with either two 
or zero electrons in the second energy state, then 
the n - 2 band would be either full or empty, 
respectively. In either case conduction would be 
impossible. Consider the case in which the n = 2 
band is empty. If the forbidden gap was made small 
enough, some of the atoms in the n = 1 band might 
gain enough energy to jump into the empty band. The 
filled band is called the valence band, and the empty 
band is called the conduction band. When the gap is 
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small enough, some conduction occurs. By 
decreasing still further the gap between valence and 
conduction bands, the conductivity is increased. 
Finally, when the bands overlap (zero forbidden gap), 
the situation is again metallic as in the case when the 
band was half-filled. Figure 3 shows the four general 
cases. 



a d 

CONDUCTOR INSULATOR 

Figure 3. Four general 

Any model chosen for calculations on conductivity 
should yield a band structure varying as a function of 
the parameters involved. This requirement rules 
out the Fermi -Dirac method, and as mentioned 
earlier, the Kronig-Penny Model is not convenient 
when lattice spacing is the parameter of interest. 

To retain the band structure, a periodic lattice must 
be used, and a sinusoidal potential with appropriate 
assumptions is chosen as the model for calculations. 
Because of the fact that symmetry is to be neglected 
and a general crystal will be used rather than an 
exact model of Mylar, a one -dimensional chair of 
atoms will be used for calculations. Such a model 
should qualitatively agree with a real crystal if the 
approximations made are reasonable. 

Figure 4 is shown to emphasize that unlike a 
single atom ( Fig. 3a), there are no free energies 
extending over all positive energy states. The free 
electrons in a crystal still have forbidden energy 
regions. This is perhaps best understood by 
remembering the wave nature of the electrons 
propagating through the periodic structure of the 


crystal. The crystal may be considered to be simply 
a filter that will pass certain frequencies. For those 
values of electron energy that are totally reflected by 
the crystal, there are forbidden gaps. Under normal 
pressure these forbidden gaps never disapjjcar 
completely, no matter how much the energy 
increases [ J . 



c d 

CONDUCTOR SEMICONDUCTOR 


types of band structure. 



Figure 4. An illustration of the difference in discrete 
and continuous energies of single "atom 11 and 
continuous energy bands and forbidden energies 
of periodic potential. 

The first potential chosen was the simplest one 
that would result in a lattice periodic in x with a 
as the period; i.e. , 

V(a) = V 0 (a) cos . ( i) 
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Since the primary interest is in the effect of 
changes in o/ t it is desirable that the average value 
of energy change as a function of a. Otherwise, no 
change in the average energy of the energy bands will 
result. This is seen by averaging the potential over 
one cycle. a 

f V 0 (a) cos — dx 

<V(a)> = -2-r =0 . 

N ' avg a 

f dx (2) 


notai, like any other solid, must be treated in 
quantum mechanics ns a gigantic molecule subject 
to the exclusion principle as well as to other 
quantum mechanical laws. 





The result is seen to be independent of a. The 
situation remains unchanged for any multiplicative 
factor not a function of position, x. In order to 
obtain a potential whose average value depends on the 
parameter o', an additive term must be used. The 
new potential may be represented by 


ION (+e) A10M ION (+e) 


Figure 5. Model chosen to determine V^cv) . 


V( a) = v,(a) +V 0 (or) cos , (3) 

with average value a 

f Vj ( a) dx 

< V < a )>avg = a = V ‘ (Q,) • 

/ dx (4) 

0 

To find a general potential,applicable to a wide 
variety of materials, an approximation must be made. 
For the present it is assumed that each atom in the 
crystal has one conduction electron. For alkali and 
monovalent metals, there is direct experimental 
evidence of this. The conduction electrons form a 
gas different from an ideal gas because of interactions 
with the lattice and with each other. These conduction 
electrons will be considered as relatively free elec- 
trons to the extent that the core can be considered 
as a separate entity with Z-i electrons and Z 
protons in the nucleus. Gauss's theorem states 
that a spherical charge distribution as seen from 
outside the sphere can be considered to be concen- 
trated at the center. This results in a charge 
+Ze + (Z-l) (-e) = -te at the center of each atom. 

The model chosen to determine V^a) will consist 
of one atom surrounded by two singly charged 
positive ions (Fig. 5) . The free electrons from the 
ions will be considered as uniformly smeared through- 
out the crystal. The method of treating electrons 
as moving freely through a periodic field has been 
criticized because it neglects collisions between 
electrons. The answer to this criticism is that a 


The next step is to determine the potential energy 
resulting from this configuration. For continuous 
charge distribution the potential energy is described 
by [6] 

V = JL ff P(x) P&) d«xd»x' . (5) 

Ix-X'l 

where p represents charge density. These 
densities are found by use of the quantum mechanical 
wave_functions. If an atom is in state n, then 
1 ( r ) | 2 d t is the probability that the valence 

electron is in the volume element d r at the point 
(r) . An equivalent statement is that -e | \j> ( 7 ) | 2 

is the average charge density resulting from the 
valence electron at the point (r) [7] . 

The actual wave functions for atoms have a 
number of nodes near the nucleus. For many 
problems, including this one, these nodes are 
unimportant and can be omitted with no loss of 
generality. Slater [8] has shown that the radial wave 
functions are given in the above approximation by 

-k r 

_ .. n n 

R (r) = r e * , 

n 

where 

k = (Z-s)/n, 
n 

Z = atomic number, 
s = Slater's Atomic Shielding Constant, 
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and 


n principal quantum number. 

The complete wave function is a product of the 
radial wave function, the angular wave function, and 
a normalizing factor, N^(r), where i and m are 

angular momentum quantum numbers. Spin is 
neglected. 

0 n U) = N n (r) R n (7) Y™ (r) . (6) 

The energy between the atom and one ion is first 
calculated. The situation is shown in Figure 6. 



Figure 6. Coordinate system for ion-atom pairs. 

For the coordinate system shown in Figure 6, 
the wave function for the ion is given by 


<l>j ( r 2> 

1 

= 6(COS ° 2) 6(<W 6(r2 ' C) Y » ( * 2 ) ’ (?) 


f <l>* <l‘ n < r 7) 

all 

space 

f 2n 2 

L'j exp { -2k^ i*j) Tj dr 1 




, A v 

Y ( r j) dfi, 


(9) 


The angular integral is equal to 6, , 6 or 

i , i m , m 

unity and the radial integral evaluated at the limits 
gives the value 


/ Ui) >P n (rt) dr t 
all 

space 


(2n+2) 


< 2k n> 


2n+3 



Thus the normalized atomic wave function is 


( 10 ) 


r" exp(-k n r 1 ) Y 1 " (rj) . 

(ID 

The values of k^ vary with Z and with the 

electron under consideration. A complete treat- 
ment of the atomic shielding coefficients is given 
by slater [8]. Values predicted from these 
coefficients generally agree with experimental 
values within 10 percent, which is close enough 
for most qualitative predictions. 

The wave functions derived satisfy the charge 
condition 


= (») 


nfl 


2k 


[2(n+l)J \ 


f pdT = -he f ip 3i< ip d t = ±e 

all 

space 


( 12 ) 


where 1 / a is the normalizing factor obtained from 


f <<r 2 ) V r 2 ) r 2 dr 2 dfi 2 = c? . (8) 

all 

space 


The interaction term l/|r 2 - r 4 | can be 
expanded in spherical harmonics as follows [9] 


+x 

l 


4 IT 

2X+ 1 




(r 2 ) 


(13) 


To obtain the normalizing factor for the atomic 
wave function, the same integral is evaluated as 
follows: 


where r < ( r^) is the lesser (greater) of r t and r 2 . 
Equation (5) becomes 
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, ,<< r ') *i < r *) V 2) 

r / — — — drj dr 2 


I i*i * i’2 I 

v v /i»V /«l_\ 
a a'» ,h '■ a \ / V aA "^ <*"* ! 


(14) 


/' r| dr, /’ r* ilr, A ( r* - ■*)r; ,n «‘X|i( -ak r,) — jy-j 

ii o r _ 

<i a 1,1 1 * m 

f f VS (r 2 ) «,) A(fc) Y 0 ( ?,) A tf,, Y x A (^, Y| m (r,) Y™!?,) .IU, Jli, . 


(15) 


for X 0 and A 1. For A 2, another possibility 
presents Itself. However, for the moment, the model 
will bo restricted to H angular momentum stntos 
( i”o) . The general expression is 


v - j? y y J5jq a _2£ti l_ 

2 VO mt'.A (2,H2) ’ (2X41)V* 


-oxp(2k n 0f) 

(2kn) A 


(2k «) 


( 2n + A + 2) 


I (2tH A ^ 2) *, 

' J (2k >* 


x C(i , i ,A; mm m. ) C(i, t 0 0 0) Y (cos tfj = 0, - 0) 

* A (18) 


The upper limit on the integration in the atomic 
wave function space was chosen to be the lattice 
parameter a. This is obviously a maximum. Later 
this may be decreased to some optimum value. The 
angular integral after integration over ft 2 becomes 


A sample term in which n=2, i=0, m=0, A=0, and 

m =0 is as follows 
A 

5 4 . „ . +720 

m=0 (19) 


n c 2 
I =0 


exp(-2k “) 
n 

2(720) 


(2k a) B + 6(2k a) 


Angular Integral 
1 m A 

— Y a (cos () 2 - 0, <p2 — 0) 


/ 


A . A . a . __m a . 

( r j) Y (r t ) Y (r t ) d 


(16) 


Expressing the integral in terms of Clebsch- 
Jordan coefficients, the complete angular integral 
becomes 


The curves for n=i, 2, and 3 are shown in 
Figure 7. In interpreting the curves, it must be 
remembered that for extremely small values of cv, 
the above model does not even approximately describe 
the situation. However, for the values of the lattice 
parameter that are of physical interest in this 
problem, it is hoped that the behavior is realistic. 
Since the situation is such that this energy opposes 
the binding energy of the electron, the sign is 
reversed when it is added in the equation. Except 


Angular Integral 
m 

= ( 4^ Y x A (cos 0 2 = 0, <p 2 = 0) 


2i + l 


[4 ir(2\+i)] 


x C(f , l , A;m, m, m ) C(f,f,A; 000 ) 
A 


(17) 


In order not to consider Clebsch -Jordan 
m A 

coefficients for which Y vanishes. Table 1 was 
A 

set up for values up to n=3. From Table 1 it is 

seen that only Y 0 , Y 2 and Y 2 contribute to the 

angular integral for values of n up to n=3. To 

evaluate the Cleb sc h -Jordan coefficients, a specific 

value of n is chosen. Choose first n=l. 

Obviously i =0 and m=0, and the Cleb sch -Jordan 

coefficient is C(00A; 00m ) C(00A; 000) from 

A 

which A=0, m^-0, and the term considered is 

unity. For n=2, the same situation applies because 
as is seen from Table 1, the values for \=1 are 
multiplied by a zero term. For n=3, the same holds 



PHYSICAL VALIDITY 

Figure 7. Additive energy as a function of 
lattice parameter. 
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TABLE 1. CLEBSCH -JORDAN COEFFICIENTS 
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±1 
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2 

l2 

2 
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4 \ 2tt / 
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for extremely small (nonphysical) values of the 
lattice parameter, this model gives exactly the 
behavior that would tend to force the n=2 band and 
the n-”3 band into conduction. The next step is to 
calculate the bands. This is done in the next section. 


The exact form of V 0 ( cv) will be discussed 
later. Now Schroedinger’s equation assumes the form 


*1 

2m 


v 2 4> +[ -V t (a) + V 0 ( a) cos 1 i|» - E ip . 
X L J (22) 


In quantum mechanical problems the energy is 
generally described by Shroedinger's equation, 

~ V 2 ip+Vip = Eip . (20) 

2 m 


This can be recognized as the Mathieu equation, 
which has been treated quite often in the literature 
[ 10] . The standard form of the Mathieu equation 
is as follows: 


Here V* is the potential and E is the energy 
associated with the particle described by the wave 
function ip. The classic case is, of course, the 
Coulomb field of the hydrogen atom. As mentioned 
earlier, however, the crystal under consideration 
can be regarded as a gigantic molecule to which the 
laws of quantum mechanics must apply. The change 
in Schroedinger's equation must occur in the potential 
V representing the periodic field of the crystal. For 
simplicity and for reasons already discussed, this 
will be considered one-dimensional. Thus 


^ a - 2 qcos 2z ^ \p = 0 ( 23) 

Letting 7tx / ot = z, equation (22) becomes 


4 [Wim _ 2mo? 2 V Q (oj_ 
[ U*) r 


) - 0 


This is in the standard form of the Mathieu equation 
when 


V = V t ( a) + V 0 ( oi) cos 


2ttx 

a 


( 21 ) 


2ma 2 [E^ y aa)) Y±i*), 

<«*.)* ’ q (TrXP 


(24) 
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Thq solutions to this eigenvalue equation yield the 
Mathieu functions of even and odd order and the 
characteristic or eigenvalues a. The Mathieu 
functions are of little interest in this problem. The 
eigenvalues a, however, yield the energy 
eigenvalues E for the electrons in the periodic 
potential. The equation at first glance appears rather 
simple, but the appearance is deceiving. Power 
series solutions yield eigenvalues in the form of 
infinite continued fractions. Various approximations 
are valid at the origin and in asymptotic regions. 

A polynomial in q was used for q very small until 
it was found that the expression was not valid over 
the range of q required by the physical problem. 
Rather than attempting to use different approximations 
in different regions, matching them where they 
overlapped, it was decided to make use of standard 
tables compiled by the National Bureau of Stand- 
ards [ 11] . Although the tables are for a slightly 
different standard form of ±e Mathieu equation, they 
are easily convertible to the form of equation (23) . 

For purposes of calculation, the one -dimensional 
crystal is considered infinite in extent. Thus, 
boundary value aspects are avoided. One would 
naturally suppose that eigenvalue solutions for a 
periodic potential neglecting boundary conditions 
would be independent of position x. This is the 
case, and it results that the nth eigenvalue is 

solely a function of q (which is in turn independent 
of x) from equation (24) . 




Figure 8. Slater representation of crystal potential. 
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Figure 9. Potential energy of electron in hydrogen 
and alkali atoms. 

It is shown from Figure 10 ( and other combi- 
nations) that to a first approximation 


= 

n 2ma* 


a + V t (o) = V^a) + C( a) a [q( or)] 


(25) 


q (at) = 


rna* 


(irSCr 


V 0 («) 


[ at | |V 0 | = constant * 35 
therefore 


(26) 


V 0 <a) » 


35 

a 


As is readily seen from equation (25) , much of the 
problem hinges upon the choice of V 0 ( a) . Slater [ 12] 
has chosen an a 2 dependence for V 0 . This is 
shown in Figure 8 when it is seen that the barrier 
height varies as a 2 . This representation make use 
of parabolic potentials that are quite often used as 
approximations in physics. At low pressures this is 
probably a good approximation; however, at high 
pressures (small o') this seems not to fit the case. 
Rather than a parabolic well, the actual case is that 
the potential goes toward minus infinity at each 
nucleus. The potential for a single nucleus is shown 
in Figure 9 [7] . By superimposing potentials of the 
kind shown in Figure 9, a new dependence on a 
arises. This is illustrated in Figure 10. 


This will be the potential used for the change in 
potential resulting from a change in lattice parameter. 

There is another factor entering into the change 
of potential with lattice parameter. This is the 
factor of ionization. As more and more energy is 
imparted to each atom (as signified by a decrease in 
ot and a corresponding rise in temperature) , the 
atoms lose their most energetic electrons and become 
ionized. The ultimate case is, of course, a plasma, 
but intermediate cases need to be taken into 
consideration. Ionization in this case must be 
treated statistically because of the complexity of the 
situation. The first factor introduced to take into 
account the ionization is the simplest possible factor 
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INDIVIDUAL POTENTIAL 
AVERAGED POTENTIAL 



Figure 10. Variation of potential as a function of lattice parameter. 


that will produce the correct action. An optimum 
value of the lattice parameter is chosen to be a 0 , 
and the ratio of at 0 to the lattice parameter is used 
as a measure of ionization. The ratio is cubed for 
more reasonable results. 


Ionization factor = 



(27) 


For example, let sulphur have a 0 = 5.0 angstroms. 
Thus for a = 5. 0 the ionization factor is unity. 

For a? w 2.0 angstroms, the ionization factor is 
approximately 16.0. For very small values of a, 
the ionization factor obviously becomes too great. 

It has already been mentioned, however, that at these 


values the additive potential is also inaccurate. Con- 
duction occurs before these values are reached, so 
the model need not be accurate to the same degree 
for very small lattice parameters. One further 
desirable feature of the above ionization factor is 
that at large values of the lattice parameter, the 
factor contributes a sort of shielding that is expected 
in such a circumstance. As the atoms separate, they 
gradually lose even their singly ionized character so 
that they would be better described as a Van der Waal’s 
gas. This situation does not arise in the problem 
under consideration, but it is helpful to know that the 
ionization factor chosen works not only in a specific 
area, but also over a range of cv's. At this stage 
then, it is seen that the potential V 0 is described by 
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V 0 ( a) 



or 


q 


rv 


35 m Qfj 


(28) 


(29) 


As can be seen in Figure 9, the width of the bands is 
related to the principal quantum number n. (It 
should be realized that this quantum number, although 
referring to levels resulting from a central potential, 
has meaning when thought of in terms of the limiting 
states in the limit in which the crystal expands into 
isolated atoms.) The widths are obtained in the 
model by using ionization factors similar to 
equation (27) with n-dependent exponents (and, of 
course, radii) . Thus, 


n-band ioniz a tiorV factor = 



(30) 


At present the exponent f(n) is determined 
empirically, although a theoretical basis is being 
considered (Fig. ii) . A further approximation is 
made for the energy bands as expressed in equation 
(24). The parameter ot is set equal to a 0 . 
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Figure il. Illustrating n-dependence of bandwidth. 

E = a - V^a) . (31) 

n 2mQ!j n 1 ' 

This does not affect the point where energy bands 
intersect. It does slightly affect the shape of the 
curves in these regions. The approximation is 


being studied in order to give a better explanation 
or a better approximation. Using for the above 
approximations f( 1) - 10, f(2) = <5, f(3) = 2, 
ot\ = 1.0, <x 2 - 2.2, cv 3 = 5, 0, and rv 0 = 5.0, the 
energy band diagram is shown in Figure 12. Tho 
band structure changes with lattice parameter so that 
the solid goes into conduction for small enough 
parameters (high pressures) . Using the scheme of 
Figure 3, the information in Figure 12 can bo 
presented in a form comparablo to snapshots of the 
energy bands for a given value of the lattice 
parameter. This is shown in Figure 13. 


LA nice f'ARAMl TIR , <. 



Figure 12. Energy band structure for model used. 

The choice of ionization factor used in the 
previous calculations was made on the basis of 
simplicity and convenience while possessing a 
reasonably accurate dependence on lattice parameter. 
An attempt will be made here to better this first 
approximation. A caloric equation of state similar 
to that of water will be used. Water, which is the 
most important medium with approximately separable 
energy has an equation of state somewhat resembling 
a perfect gas [ 13] 



where A, B, and y are practically independent of 
entropy. This fact assures that by neglecting 
entropy changes, no great error will be introduced 
into the ionization factor. For reference purposes 
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Figure 13. Alternate method of depicting energy gaps resulting from model. 


Po - density of H 2 0 at 0° C , 
A -= 3001 atm , 

B - 3000 atm , 

and 


Y 


7 


A change of variable from density to lattice 
parameter yields 


P(<*) 



- B 


(33) 


At constant entropy, the work done by pressure, 
p, on volume, V, is pV; therefore, 


dW _ pdV Vdp 
do do do 

The atomic volume Is set equal to 


(34) 


The change in energy as a function of lattice 
parameter becomes 


■dW = (A-y) 



o 2 do - B o 2 do 



(37) 


However, it will be noted that to obtain one 
atmosphere of pressure under normal conditions 
(V = V 0 ) , it is necessary that A = B + 1. If the 
change of notation A = P fl is made, then the 
equation of state becomes 



(38) 


dp 


3y ( P 0 + 1) 



(39) 



and 



(35) 


(36) 


The expression for the energy change becomes 
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where the fact that da is negative has been taken 
into account. This elation can be rewritten 


dW - P 0 




(40) 


It will be a good approximation to consider 
i/p Q s o and use 



as the energy to be used in further calculations. 
Having found this energy change as a function of 
lattice parameter, it is reasonable to assume that 
as the energy change is greater, the probability is 
higher that an atom will be ionized. The simplest 
type of dependence given by assuming that the 
probability of energy input W ionizing the nth band 
electron is 


This probability can easily be integrated to obtain 


IMtt) 


f ■*. r 


0 ' fe) r (V) r 


j» 

ii 


('») n* 


(y - 0 


«y 




and, rearranging terms. 



P n ( a) = n 2 E 0 -1 -f P 0 a 3 



(46) 


This quantity is negative; however, by taking the 
absolute value of W it is seen that a reversal of 
sign occurs such that 


P n ( a) = n 2 E 0 _1 j- P 0 a 3 



(47) 


This function yields zero probability of ionization for 
zero change in lattice parameter (no compression) 
and increasing probability as the material is com- 
pressed. The probability is chosen to be maximum 
when the input energy is equal to the ionization 
potential of the nth band electron. The corresponding 
value of or is found by solving 


W ( a) = E 

n 


_ JWi _ JW 


En 1 = n ’ W » 


(42) 


or 


3y _ 3y 


= a 


1 + 




(48) 


where 
W 0 = 

and 

W= f dW ( a) . (43) 

“o 

It can be seen that the probability of ionizing 
the nth band electron has the form 

a 

P (a) =+n 2 E 0 -1 f f(a') da' , (44) 

n J 

or specifically 


jm 

|E 0 | 


(45) 


where 



which must be solved by graphing or iteration for 
each particular initial value The value of a 

found by the above procedure is the lower limit 

beyond which the probability becomes meaningless 

[P (a)> i] . The expression P ( a) should replace 
n n 

the ionization factor (a^/a) F(n) given in 

equation (30) . Results of calculations with this 

ionization factor arc not given because of a shortage 

of time and the lack of a plotting machine. It is to 

be expected that some modification of this factor 

may be necessary to prevent the ionization potential 

from disappearing when a = a . However, such a 

o 

vanishing factor is more appropriate to a dielectric 

at or = Of than a metal. It should be mentioned 
o 

that while it is conceptually useful to regard each 
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band as deriving from a particular atomic state (in 
the sense of the tight binding method) , it is only 
rigorously true if the energy band is narrow and 
separated from neighboring bands by large energy 
gaps. If several bands overlap, they must all be 
regarded as deriving from the same group of atomic 
states. This adiabatic device to label crystalline 
bands in terms of bound atomic levels obviously 
fails for all higher bands that arise from the 
continuum. In three dimensions the order of the 
levels in energy can change as the lattice parameter 
is varied. The device may then fail even for low- 
lying bands [5] . Nevertheless, the use of this 
concept to derive an ionization factor is more 
satisfying intuitively than simply choosing one which 
yields the proper results. 

Before entirely dismissing the results of the 
first ionization factor, equation (30), it should be 
noted that one of the empirical features required for 
correct results has a basis in theory. The device of 
using discrete initial values ai^ for the various bands 

can be justified by the following arguments. As 

mentioned in the description of how energy bands 

arise, the interaction of electrons initially in the 

same state results in each of the initial states being 

slightly displaced so that the final configuration 

consists of two (or more) new states. The phrase 

'’interaction of electrons initially in the same state" 

assumes an importance not immediately obvious. 

Consider two atoms separated by a distance small 

enough that the outer orbits overlap and also the outer 

orbits partially overlap the inner orbits as in 

Figure 14b. The electrons in the outer orbits and 

those in the inner orbit arc not "initially in the 

same state." Therefore, although there is an 

interaction shifting the energy levels, the levels do 

not necessarily spread out to form a band. Only 

when the separation of atoms is such that the inner 

levels overlap (ff<2xr, . ,) do these levels 

inner level 

spread into a band. This is the justification for the 

use of the a and a possible need for similar 
n 

provision in P ( cv) such that P (Qf)-*P (a ). 

n n n n 

The problem was also considered in general 
thermodynamic terms in an attempt to impose 
limitations upon the usefulness of the model. The 
treatment follows Seitz and Turnbull [ 14] and will 
not be reproduced in this review. The results can 
be summarized by stating that the change in energy 
is found to be proportional to the change in volume; 
i.e., 





<c) 


u<a m <a ° 

(d) 


Figure 14. Different possibilities of overlapping 
wave functions. 



This relation can be used to obtain an order of 
magnitude approximation to the energy changes 
involved for specific pressures. When the 
equipartition theorem is used, it follows that about 
one-half of this energy should contribute to the 
potential energy. The volume can be expressed in 
terms of the lattice parameter a (assuming an atomic 
radius of oi/ 2) as 


Vn =— tt 


cr 

2 


( 50 ) 


For a pressure of 2 Mb, the lattice parameter is 
halved; therefore 


Vi 


1 / \ 
2 \ 8 / 


and 


AV = v t - v 0 - ^ a 3 


The potential energy change is seen to be 


AE AE=-fp AV — ~ Mb 

pot 2 2 av 2 


7 o? 
16 


"TT (51) 


AE = p AV 
*av 


(49) or approximately 
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AE , - i Mb X 
pot 


atomic volume 
16 


(52) 


The pressure in megabars can be converted to 
energy per atomic volume, and the result is found 
to be 


1 Mb « o.8 


eV 

atomic volume 


which yields a change in potential energy, 


(53) 


however, mean that the pressures cannot bring about 
substantial increases in conductivity. It means only 
that a simple, linear treatment of the shock waves 
is insufficient. Interaction and reflection of nonlinear 
waves lead to enormous increases in pressures as 
has been verified experimentally. Such nonlinear 
interactions could veiy easily result in conduction 
band energy changes. They also make likely the 
appearance of hot spots wherein the required 
energies are available. The study of nonlinear proc- 
esses seems necessary to gain much more informa- 
tion about conduction brought on by high pressures. 


AE “ °‘; ' eV “ 0.05 cV 
pot 16 


(54) 


CONCLUSIONS 


Thus it would appear that at the pressures present 
in the material, there would be insufficient change in 
potential energies to affect the conduction bands 
even assuming all of the energy went into this 
direction. 

One way in which the required energies may be 
achieved is to use the super-position properties of 
interfering waves. It is quickly seen that even the 
added pressures do not yield sufficient energy changes 
to affect conditions appreciably. This does not, 


Qualitatively, the model results in conductions 
with increasing pressure expressed by means of 
decreasing lattice parameter. The thermodynamic 
analysis indicates, however, that a linear treatment 
of the problem is not capable of determining quan- 
titative behavior with acceptable accuracies. If 
carried several steps further by a statistical 
mechanical treatment, the relaxation times for the 
conducting dielectric could be computed. This 
information may be sufficient to determine meteor 
impact velocity from empirical data. 
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HIGH TEMPERATURE EQUATION OF STATE 

By 

R. J. Nau inarm 


FOR ALUMINUM 1 


SUMMARY 

An equation of state is developed that is capable 
of describing metallic elements in both solid and 
liquid phase from ambient temperature and pressures 
to states of extreme temperature and pressure. A 
novel technique for treating the atomic vibrational 
contributions is introduced by assuming each atom 
vibrates independently in a tan 2 potential. The 
Schroedinger equation can be solved exactly for this 
potential function, and the quantum mechanical parti- 
tion function is computed directly. In this manner 
a continuous, thermodynamically consistent descrip- 
tion of metals is obtained which reduces to an 
Einstein solid at moderate temperatures and to an 
ideal gas in the limits of high temperatures and/or 
large volumes. 

Detailed numerical computations are carried out 
for aluminum. All the constants required to 
configure the model arc obtainable from elementary 
thermodynamic data at ambient conditions. Excellent 
agreement is obtained with experimental shock 
compression data on solid and porous samples at 
pressures to 5 Mb and temperatures to 25 000°K. 

INTRODUCTION 


Substantial progress has been made in recent 
years in understanding the behavior of metals at 
extreme temperatures and pressures. Shock com- 
pression techniques have achieved pressure measure- 
ments to approximately 10 Mb and temperatures of 
tens of thousands degrees [ 1] . At the lower pres- 
sures, the Debye or Einstein model for energy and 
the Mie-Gruneisen equation for pressure is adequate 
and zero-degree isotherms were extracted from shock 
measurements by using these relations to subtract out 


the thermal contributions. It was found that a Morse 
potential predicted the zero-degree isotherms to a 
fair degree of accuracy [2], which confirmed an 
earlier suggestion by Slater [3] . 

As experimental pressures and temperatures 
increased, various inadequacies in the theory became 
evident. Electronic contributions were introduced. 

For compressible metals, anharmonic terms became 
significant. These were treated by Al’tshuler et al. 
using the free volume theory of Lennard-Jones and 
Devonshire [4]. Pastine [5] used the perturbation 
method of Liebfried and Ludwig [G] to correct for 
anharmonic effects. 

For the high-temperature states reached by 
shocking porous samples, it was found that the 
behavior became more ideal gas-like. Kormcr [7] 
proposes a set of empirical interpolation equations 
to transform solid-like behavior of energy, pressure, 
and heat capacity to the ideal gas relations. Urlin [8] 
proposed an empirical free energy function to 
account for phase transitions. All of these empirical 
relations require adjustable constants which are 
evaluated from experimental high-pressure data. 
Furthermore, there is no guarantee of thermody- 
namic consistency between them 2 . 

An equation of state for metals is developed in 
this work that is capable of describing the liquid-dense 
vapor phase as well as the solid phase from ambient 
conditions to temperatures and pressures exceeding 
the present experimental range. Rather than use in- 
terpolation equations to transform the behavior of a 
solid to that of a gas at high temperatures or at low 
densities, the approach will be to start with empirical 
interatomic potential functions and develop the entire 
equation of state using quantum statistical mechanics. 
The only empirical constants required are obtainable 
from elementary thermodynamic quantities such as 


1. This report was based on work submitted in partial fulfillment of the requirements for the degree of 
Doctor of Philosophy in the Department of Physics in the Graduate School of the University of Alabama. 

2. Kormer's pressure, energy, and heat capacity were chosen in such a manner to be consistent among 
themselves. 
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heat of vaporization, compressibility, entropy, etc. 

In this manner the behavior of a metal throughout the 
entire high-pressure, .high- temperature regime can be 
predicted from its properties at ambient conditions. 
Since all the thermodynamic properties are derived 
from the partial function, thermodynamic consistency 
is guaranteed. This allows a complete thermodynamic 
description , so that phase transitions are described 
naturally from Gibbs energy considerations. 


relation 

p= - OF k /av) T 


and is 



X(e 2b,1 - JI -. b,1 - £) ) 


( 2 ) 


FORMULATION 


Vibrational Terms 


It will be assumed that the Helmholtz energy can 
be expressed as the sum of three independent 
contributions: the lattice term, F^, the ionic 

vibration term, F , and the free electronic term, F . 

v e 

Since. ;all thermodynamic functions are derivable 

from the Helmholtz energy, it immediately follows 

that they all can be expressed as sums of these three 

contributions. The fact that there are interactions 

between the free electrons and the lattice will be 

accounted for by introducing an effective electron 

mass. 


At moderate temperatures, the ionic vibrational 
contributions may be computed from the assumption 
that the ions behave as harmonic oscillators. This 
assumption leads to the Einstein model, in which it is 
assumed that all atoms vibrate independently at the 
same frequency, or to the Debye model which considers 
a distribution of normal mode frequencies. At very 
high temperatures, the vibrational amplitudes are 
such that nonlinear restoring forces must be con- 
sidered. This is accounted for by assuming each atom 
vibrates independently 3 in a potential well given by 


2a 2 mcu 2 E 

*(x) =— ^ 



(3) 


Lattice Terms 


Since the Morse potential has been found to be a 
reasonable empirical representation for the lattice 
energy of a metallic solid, the Helmholtz energy 
for the zero-degree lattice is written 



where L is the zero-degree heat of vaporization, 

0 l/3 

or total binding energy, £ = (V/V ) 1 and b is an 

empirical constant determined from the compres- 
sibility. From the relation F - E - TS, the internal 
energy, E^, of the zero-degree lattice is identical to 

F^. The zero-degree isotherm is found from the 


When the displacement x is small compared to the 
atomic spacing a, the potential reduces to 

niWg 

*(x) = — 1 x 2 , (4) 

where a; is the Einstein frequency. For a harmonic 

oscillator, o> 2 = K/m, therefore, equation (4) has the 
form l/2 Kx 2 where K is the spring constant 
corresponding to the Einstein frequency. 

At larger displacements, x~±a, the potential 
approaches infinity which describes a hard sphere 
collision between point masses. This roughly cor- 
responds to nuclear collisions between neighboring 
atoms whose position expectation values are ±a. 

The choice of the tan 2 functional representation 
was made because it behaves as desired in the limits 


3. The assumption of independent vibration is justified because even at moderate temperatures the Einstein 
and Debye models give almost identical results. 
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and also because it allows a closed -form -eigenvalue 
solution to the Schroedinger equation. The energy 
eigenvalues are 4 


e n ^ n2 h € o + + ^o )^ 2 9 

(5) 

n « 0, 1,2, . . . 

where is the degeneracy energy A 2 # 2 / (8 m a 2 ) . 

The vibrational partition function Z for N 
atoms, each with three degrees of freedom, is 


In Z = 3N In £ 
n=0 


-e /kT 
e n 


( 6 ) 


where k is the Boltzmann constant. 

The various thermodynamic functions are 
obtained from the partition function in the usual 
manner. 


« _ e /kT 

F = -kT In Z = -3NkT In / e n , (7) 


n=0 


E -• kT 
v 


8 In Z 
8 T 

3N 

f,E v\ 

3N 

i,T /v “ 

T?P 


( 8 ) 


<'»> ’ 


and 


P = -kT In Z = 3N / P \ , 

v c)V \ n / 

where 
P 


8e 


av 


( 9 ) 

(10) 

(ID 


Differentiating equation ( 5) , the P contribution 
becomes 




y*V+-f- c 2 

o O 

(KW+e 2 ) l/2 


( 12 ) 


where y is the Gruneisen ratio, -(3 In u/d In V) . 
Formulation of y for solids has received extensive 
treatment in the literature. 


The brackets ^ ^ donote the ensemble average; eg. , 


< X ) "■* 


V -c /kT 
2 x o n 


n-() 


« -c /kT 

r n 
e 

n=0 


(13) 


At solid densities, Mco » e and the Boltzmann 

o 

factor reduces to exp ^-(n 2 e Q + n)fco)/kTj . Unless 

kT » Xu, the Boltzmann factor will become negligibly 
small at low enough values of n to prevent the n 2 € Q 

term from contributing significantly. In this case, 
the Helmholtz energy reduces to 


F = - 3NkT 2 In f : p _nJ[< ' w/kT 
' n=0 

= 3NkT In (i-e _ ^ W ^ kr ) , 


(14) 


which is the well-known result from the Einstein 
model . 


The Liquid Potential Function 


The zero-degree potential energy in the Liquid 
phase is represented by 


Li e 


F k = 


BV 

2b( i-|) _ o 
V 


I -£l 


(15) 


( e 2 b (l-« / S,)_ 2<5 b(l- S / f,i) tE ;{S5i 


where Lj, B, E^, and are empirical constants 

peculiar to the liquid phase, and and b are the 

same as in the solid potential, equation ( 1) . This 
particular form was chosen for the following reasons. 


4. See problem 12 in D. ter Haar, Problems in Quantum Mechanics, Academic Press, New York, 19C0. 
Actually, there should be a l/2 added to n in both terms corresponding to the ground state. However, 
the ground state energy is included in the zero-degree lattice terms; therefore, it is suppressed in the 
vibrational terms. 
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For £ — the repulsive term is the same as in the 

Morse potential, but the attractive term is the form 

resulting from van der Waals forces. This potential 

function cannot be extended to small £ because 

eventually the attractive term will override the 

repulsive term. To avoid this, tho form given for 

£ was chosen. This is identical to the Morse 

curve for the solid except that E and have 

m 

been introduced to. account for the heat of fusion and 
different compressibility of the liquid phase. 


definite position in the lattice. In a liquid, particles 
can exchange positions and thereby lose their 
distinguishability. The partition function must be 
adjusted accordingly for ’’proper Boltzmann counting.” 
Since each atom occupies a coll with volume (2a) [3], 
n volume containing N atoms spaced at average 
distance, a, can have N/s colls. There are (N/8) 
ways of arranging N particles among N/8 cells, of 
which N! arc redundant because of particle indistin- 
guishability. The partition function becomes 


The pressure is again given by = - ( dF^/clV) 
and is 


p k = 


2bU_ 2b( i-|) 

3 t 2 y e ’ * *1 

s o 


2bL 0 

3£l? Z V 




( 10 ) 

: • 


Of the four additional constants introduced in this 

potential, two are required to match the F, and P, 

k k 

at £ = {j. The remaining two are chosen to produce 
the observed energy and pressure at the melting 
point 5 . Figure 1 compares the liquid potential with 
the Morse potential for aluminum. 


Z - 


(N/8) N 

n: 


3N 


oO 


? 

n=0 


-e /kT 
n 
e 


(17) 


The Helmholtz energy becomes, using Stirling’s 
approximation. 


oo /kT 

F = NUT (In 8 - 1) -3NkT In J e n 
v n=0 


(18) 


The energy eigenvalues, are given by 

equation (5) ; the internal energy and heat capacity 
are given by equations (8) and (9), respectively, and 
the pressure is given by equations ( 10) and ( 12) , 
except that the Einstein frequency and the Gruneiscn 
ratio will have different values for the liquid. 


The Liquid Vibrational Component 


The Lennard-Jones Devonshire model of liquid 
and dense vapors treats each atom as though it were 
moving In a cage of its nearest neighbors which are 
considered fixed. Since the number of nearest 
neighbors is usually 12, a high degree of symmetry 
exists, and a spherically symmetric potential maybe 
assumed. Following these concepts, the same form 
of potential function as assumed for the solid, 
equation (3) , is used to describe the behavior of an 
atom in such a cage. 

The partition function will be somewhat different 
from equation (6) , however, because of particle 
indistinguish ability. In the solid phase, particles 
are, in a sense, distinguishable because of their 


It is interesting to examine the behavior of the 
liquid in the limit of large volumes where Xw— -0. 

In this limit, the atom moves freely within its "cage." 
Encounters with the wall correspond to elastic point 
collisions characteristic to an ideal gas. Since 
kT » e , the summations may be replaced with 

integrations and 

OO 

l e 

n=0 


(n 2 + 2n)e 

c 

kT 


\ 1/2 
rr kT \ 


4e 


o 


(19) 


The Helmholtz energy becomes 


2_ 

Mb-T 1 n 

1 

H 

/ \ 2/3 

V_) 

3 

in k 1 in 

2 7 rtf 2 

\ N / 

- 


. ( 20 ) 


5. The parameter, , is the value of £ corresponding to minimum potential in the liquid phase. It will 
generally be more than unity because the distribution of the atoms in liquid metals is random instead of the 
more efficient close packing usually associated with crystal lattices. This choice of functional dependence 
was made to assure P (£j) = 0 and to make the F and P approach the same values as the solid phase 
at small k k 


82 



R, J. NAUMANN 



0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9 2.1 2.3 


1 = (V/Vo) 1 * / 3 


Figure 1. Comparison of solid and liquid potentials. (The liquid potential approaches the solid at high 
compressions, differs by approximately the heat of fusion at V = V , and exhibits a van der Waals 

behavior at large volumes.) 


Similarly, the other equations reduce to the ideal gas 
2 NkT 

relations, E = — NkT, P = ■■ , and 

v 3 v V 

3 

= — Nk. Using the result of equation (20) , the 
relations S= (E-F)/T, and the fact that E^= F^, 
the entropy in this limit becomes 



which is the well-known Sakur -Tetrode equation for 
the entropy of an ideal gas. 

At liquid densities, and the treatment 

of the vibrational component of liquids at moderate 
temperature reduces to the Einstein model. However, 
the Einstein frequency and the Gruneisen ratio will be 
different in the case of a liquid. 
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A crude estimate of tho vibrational frequency can 
bo made by considering on atom at tho origin with 
nn atom at ±a. The change In potential resulting 
from a displacement x Is 


grand partition function, one oblnlnH [0] 

N r e 1//z dc 

V ‘ 8 j . 


*(x) = 4>(a+x) + *(a-x) - 2*(a) (22) 

Using only the repulsive term in the Morse Potential, 
2b(i-r/a ) 

*(r) = e 0 , 

o 


*(x) = 2 * Q e 


2b( l-{) 



(23) 


and 


E = g 


r e 3/ ^ 2 de 
{ e'^ 7ifr 


where 

2 7 / 2 7rm 3 ^ 2 

g = 

h 3 


(27) 


(28) 


For small displacements, 


$(x) w 2 i e 
o 


2b( l-£) 4bV 


(24) 


This has the form of a harmonic oscillator potential 
with a frequency 


4 $ b 2 

e 2b(1 -^ 

m a* 


(25) 


from which y = — b £. The w as a function of 

u 

volume can be found from equation (25) and the value 
of at 5=1. This can be obtained from entropy 

measurements. In the limit of the Einstein model, 
the entropy becomes 


The chemical potential /z must be found for a 
given volume and temperature by a numerical iterative 
solution of equation (27) , then equation (28) may be 
integrated to find E. The pressure is obtained from 
2 

the identity, PV = — E, which holds for Fermi 

O 

as well as Bose gases. Since, by definition, the 
chemical potential is the Gibbs energy per atom, 

G = Ity. Using the identity G = F + PV, the electronic 
contribution to the Helmholtz energy may be found. 

For temperatures that are small compared to 
the Fermi temperature, the integral in equation (27) 
may be approximated 




ttVt 2 

12 


(29) 


/ XWkT \ , 

\ 7°^)' 


S = 3 Nk 


-Nk (In 8-1) 


, , , -Kw/kT 

3Nk In ( 1-e ) 


and equation (28) becomes 

/ N k 2 T 2 


(26) 


E =T Ne F + 4 


(30) 


Given S and T, the above expression can be solved 
for w. 

Figures 2 and 3 show how the vibrational com- 
ponent of heat capacity approaches ideal gas -like 
behavior in the limitof high temperatures or low 
densities. 

Electronic Contributions 


where the Fermi energy e is 

r 



The first terms in equations (29) and (30) are the 
zero-degree degeneracy energy. This is already 
contained in the lattice terms and will therefore 
be suppressed. The electronic contribution to the 
Helmholtz energy is 


The electronic terms are obtained by treating 
the free electrons as an ideal Fermi gas. From the 


i? N k 2 T 2 
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Figure 2. Behavior of heat capacity for solids and liquids in the high temperature regime. (Breaks 

in the curves represent phase changes.) 


The Fermi energy is computed from the observed COMPUTATIONAL RESULTS 

electronic heat capacity coefficient which can be 
measured at low temperatures. The C is obtained 

by differentiating equation (30) , giving 


Specification of Constants 


c 

ve 



7T 2 Nk 2 T 
2e_, 


(32) 


At very low temperatures, the electronic terms 
dominate the heat capacity, which is experimentally 
observed to depend linearly on T. The value of 

obtained in this manner will be somewhat different 
F 

from the value obtained from equation (31) if the 
free electron mass is used for m. This difference 
arises from the interaction of the electrons with the 
lattice. The effective electron mass is defined as that 
mass which makes equation (31) agree with the Fermi 
energy obtained from the experimental heat capacity. 


The equation of state developed in the preceding 
section was applied to aluminum. The various 
constants required are obtained in the following 
manner. 

The zero-degree heat of vaporization is 

found from tables [10] to be 86 400 kcal/mole. The 
volume at 0°K is obtained by extrapolation of thermal 
expansion data and is taken to be 2. 736 gm/cm 3 . 
From compressibility data, b is found to be 3.241. 

The Slater relation [ 11] , 

2 1 

y - - — - — VPy/p 1 , where primes denote 
u 2 k k 
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Figure 3. Transition of heat capacity from liquid to gaseous behavior. (The curves indicate the behavior 
of the vibrational component of the liquid phase only, and ignore the liquid -gas coexistence region.) 


differentiation with respect to V, gives y to within 
a few percent of the experimental value at normal 
density. Because of this agreement the Slater relation 
was used for the calculation for solid aluminum 
despite the objection that it assumes a volume 
independent Poisson ratio. In the development of the 
Slater relation, the result 


aluminum, the Einstein temperature, 0 - hWk was 
found to be 269. 38° K. For the liquid phase, the 
entropy in equation (26) is used to find the liquid 
Einstein temperature which is corrected to V q 

using equation (25) . For liquid aluminum, the 
Einstein temperature was found to be 117.935°K. 


w 


2 

E 




is obtained. Given to , the Einstein frequency at 
o 

V , the to at any other V can be found. The 
o E 

quantity to^ is found by comparing the entropy of an 

ideal solid in the Einstein approximation to the 

measured value at low temperature and using the 

Slater relation to correct to V . For solid 

o 


The electronic terms for aluminum are obtained 
by considering all three valence electrons as free and 
by taking the effective mass to be 1.6 times the normal 
electron mass. This brings the calculated electronic 
heat capacity, equation (32), into agreement with the 
measured value [ 12] . 

The constants L^, B, E^, in the liquid 

potential are found in the following manner. At the 

melting volume ratio £ and temperature T , the 
m m 
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energy Is E y «„, y - (y, Tj 

winch must equnl to the energy of the zero-degree 
crysLnl U phis the hent to melt, which is known 


ojt|KM*linc.n(.:illy . All quanlllUiH required to cnnqiulo 

I'M/ » 'I 1 ) and K (/ , T ) h live boon specified; 

v ‘m in e 'm m 

thus, Iv (/ ) is known. Similarly, 

k in 


\\ (« ) I I* (£ , T ) + P (£ , 
k s m v in m o s m 


T ) 
m 


atm ’ 


therefore, I* <€ ) is known. Using equations ( 15) 

and ( 1 6) for £ ^ £ ls simultaneous equations 
containing and B are obtained. Their solution 
for aluminum yields Lj = 72 53 1 cal/mole, and 
B -■ 140 397 cal/mole. The quantity £j is found by 
requiring the two expressions for P , equation (16), 


to be equal at £ = For aluminum, it was found 
that 4!= 1.0265. Finally, the requirement that both 
expressions for E^, equation ( 18) , agree at 

£ = £i is used to obtain E = 2654 cal/mole. 

* m 


The Fusion Curve and Equation of 
State Surfaces 


(laving specified the various constants for 
aluminum, computations were carried out for both 
the solid and liquid phase for various values of £ 
and T. Isotherms of a Gibbs energy versus pressure 
plot are shown in Figure 4. Since the phase that 
produces the lowest Gibbs energy at a given temper- 
ature and pressure is the stable phase, the inter- 
sections of the solid-liquid isotherms represent 
points on the fusion curve (Fig. 5) . Using the fusion 
curve, isotherms are constructed on a F-V surface 
(Fig. 6), an E-V surface (Fig. 7), and an S-V 
surface (Fig. 8) . Isoenergy lines on a P-V surface 
are shown in Figure 9. 


Comparison with Experimental Results 


To compare the equation of state developed 
in this paper with data obtained from shock com- 
pression experiments, the Hugoniot requirement, 



Figure 4. Gibbs energy isotherms. (Circles 
represent the intersection of the liquid and solid 
isotherms and determine the melting pressure for 
that particular temperature.) 



E - E = 4- P(V - V) 
o 2 oo 


(33) 


Figure 5. Melting curve deduced from Gibbs 
energy plot. 


is solved simultaneously with the equation of state. 
In experiments with porous samples, the volume 


V qq represents the actual specific volume of the 
sample, which may be several times the volume of a 
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v/v 0 

Figure 6. Pressure isotherms. (Breaks in the isotherms represent transition from solid 

to liquid phase. ) 


normal crystal. The states that can be reached 

by shocking a crystal with initial specific volume 

ratio V must lie along the intersection of 
oo 

equation (33) and the equation of state surface in 
P-V-E space (Fig. 9) . By varying the porosity of 
the sample, any state can, in principle, be reached 
by shock compression. 

In performing shock compression experiments, 
the observables are the shock velocity D and the 
material velocity u. These are related to the 
thermodynamical quantities E, P, and V through 
the Hugoniot relations 

(D-u) = VD (conservation of mass) (34) 

PV = Du (conservation of momentum) . (35) 

These relations are used to relate the observables 
u and D to the Hugoniot states found by solving 
equation (33) with the equation of state. Figure 10 


compares the predicted results for various porosities 
with measured results summarized in Table i. Very 
good agreement is obtained considering that the 
equation of state constants are determined solely from 
ambient properties and contain no constants adjusted 
to fit the high-pressure data. The highest pressure 
datum point falls somewhat above the predicted curve 
for liquid aluminum, but is below the curve for solid 
aluminum. The fusion curve clearly indicates that 
melting should have occurred at this state, but the fact 
that the material was shock compressed may have 
prevented the phase transition to become complete in 
the short time involved, and the discrepancy could be 
that the aluminum is behaving as a super-heated 
solid. On the other hand, the inaccuracies in the 
model, particularly in the assumption of a Morse 
potential configured from the ambient compressibility, 
could easily account for such a discrepancy at 
pressures of 5 Mb. 

The reproduction of states on the Hugoniot for 
normal crystal density is not a particularly crucial 
test for an equation of state. The lattice terms are 
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* (k*/t.c) 


Figure 10 . D-u plots comparing computed results 
with experimental data points. (The short dashed 

region on the V /V = 1 curve represents the 
oo o 

transition region from solid to liquid phase. The 
long dashed curve represents a superheated solid. ) 

the major contributors at all but the highest com- 
pressions, and the fact that the restoring forces 
"stiffen up" at high compressions extends the validity 
of the assumption that the atoms vibrate as harmonic 
oscillators to high temperatures. Therefore, the 
Einstein or Debye Model together with the Mie- 
Gruneisen equation and a properly adjusted zero- 
degree isotherm will give a good representation of 
the normal Hugoniot to several megabars. Addition 
of electronic terms will extend such a model even 
further. However, such models fail for the lower 
compression, higher temperature states obtained by 
shock compressing porous samples. The equation of 
state developed in this work successfully predicts 
these states, as shown in Figure 10. 

It is useful for shock compression work to plot 
various thermodynamic coordinates against material 
velocity u as a means of expressing the thermo- 
dynamic state attained in a shock process. The 
material velocity is particularly useful because it can 


bu easily related to tho relative velocity of the 
impacting samples. Figure 11 is a plot of pressure 
for dif for out porosities vorsus u, and experimental 
points are shown for comparison. 

Figure 12 is a plot of temperature in the shocked 
region as a function of u. The various datum points 
are illustrated not as a comparison, since measure- 
ment of temperatures were not made, but as an 
indication of the range of temperatures accessible to 
this technique. The fact that higher temperatures 
are produced by increasing the porosity is dramat- 
ically illustrated by comparing the temperature for 
the 1.43, 2.08, and 2.98 points and 1.00 point at 
5.62 km/sec. These states were created by the same 
impact velocity. The fact that there is a discontinuity 



U ( km/ft«C ) 


Figure 11. P-u plots comparing computed results 
with experimental data points. (No abrupt change is 
apparent that would indicate a phase transition. ) 
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TABLE 1. SUMMARY OF IIUGONIOT DATA FOll ALUMINUM 


GROUP I a 

D (km/sec) 

(km/ sec) 

P(kb) 

v/v 

O 

7.531 

3.230 

335.8 

0.7874 

6.927 

2.319 

222.7 

0.8333 

6.500 

1.700 

153.5 

0.8696 

GROUP II b 

D (km/ sec) 

■ 

u (km/sec) 

P(kb) 

V /V 
o 

9. 13 

2.80 

693 

1.442 

10.39 

3.70 

1042 ; 

1.553 

12.94 

5.62 

1971 

1.767 


GROUP III C 


Voo/^oo 

D (km/sec) 

P(Mb) 

V /V 

O 

d 

u 

1. 00 

18.31 ± . 16 

4.93 

2. 185 

9.93 

1.43 

11.74 ± . 10 

1.391 

1.498 

G. 25 

2.08 

11.42 ± .09 

1.003 

1. 176 

G. 74 

2.98 

10.75 i- .08 

0.702 

1 . 0 15 

7. 18 


JJ 

Those are representative points of those listed by J. M. Walsh, 

M. II. Rice, R. G. McQueen, and F. L. Ynrgcr, Phys. Rev. ; 108 , 

196 ( 1957) The Material velocity u is taken to be 1/2 the measured 
free surface velocity, u^. 

V. Al'tshuler, S. B. Kormer, A. A. Bakanova, and R. F. Trunin; 
Soviet Physics - JETP, li, 573 ( 1962) 

C S. B. Kormer, A. I. Funtikov, V. D. Urlin, and A. N. Kiksnikova, 
Soviet Physics - JETP, JL5, 477 ( 1962) 

^The quantity u was not stated in the paper. It was recovered by the 
relation u - [1- (V/V q ) <V 300 /V oo ) ] D. 
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u ( km/ *«c ) 


Figure 12. T-u plots. (The data points are shown 
only to indicate range of temperatures reached by 
present experimental techniques. The break in the 
curves indicates that temperature measurement 
could serve as a method of detecting 
phase transitions. ) 

in the temperature curve at the melting region gives 
a possibility of experimentally detecting the melting 
point at very high pressures, although such temper- 
ature measurement would be extremely difficult to 
make in solids. 

Figure 13 is a plot of entropy as a function of u. 
As before, the indicated data points are intended only 
to indicate range of experimental states, not measure- 
ments. For comparison, the temperature and phase 
of aluminum at ambient pressure for various values 
of entropy is indicated. Assuming the release is 
adiabatic and isentropic, an estimate of the release 
temperature can be made. Attempts have been made 
to measure the release temperature of shock 
compressed solids. Taylor [ 13] reports a favorable 



u ( km/sec ) 

Figure 13. S-u plots. (The data points indicate 
measurement range rather than actual 
measurements. The dashed lines indicate entropy 
values corresponding to the stated temperatures at 
ambient pressure. If the release is assumed to be 
isentropic, these can be used to indicate release 
temperature. ) 

comparison between calculated and measured release 
temperatures for Cu up to the melting point. 

CONCLUSIONS 

A complete equation of state for metals that 
describes both the solid and liquid-dense vapor phase 
has been developed. All empirical constants required 
can be obtained from elementary thermodynamic data 
at ambient conditions. Since all thermodynamic 
functions are derived from the partition function, 
thermodynamic consistency is guaranteed. Phase 
transitions between the solid and liquid come about 
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naturally by using the Gibbs energy to determine the 
stable phase for a given state. Detailed computations 
were carried out for aluminum and the ability for the 
model to predict the behavior of aluminum for 
pressures of 5 Mb and temperatures of 20 000°K was 
demonstrated. Since only elementary thermodynamic 
data at ambient conditions are required to configure 
the model, extension to other metals is straight- 
forward by insertion of their appropriate constants. 


The model can be extended to higher temperatures 
with some additional computational effort to obtain 
the electron contributions for temperatures com- 
parable to the Fermi temperature. The model can 
also be extended to lower densities, provided the 
temperatures arc low enough to completely neglect 
electronic contributions. 
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LOW LIGHT LEVEL TV OBSERVATIONS 

By 

Stuart C lilton 


A numlxjr of ways have boon employed to inves- 
tigate the mass-flux characteristics of the near- 
earth meteoroid environment. Visual and photograph- 
ic means have been used to define such charac- 
teristics in Lhe larger mass region to as small as 
10" 1 gm , while satellites , notably the Pegasus 
spacecraft, have extended observations of small 
mass meteoroids to as large as 10 -6 gm. The mass 
region between 10 _1 gm and 10 -6 gm is a most 
important region in regard to long term spacecraft 
such as the S-IVB workshop. This region has been 
investigated primarily with radar observations. 
However, problems and uncertainties inherent in 
the analysis of the radar data have ted to the desira- 
bility of using alternate methods to examine this 
mass region. 

Such a method is provided through the use of 
low light level television systems such as the image 
orthicon, the SBC vidicon, or the image isocon. 

For this reason, the Meteoroid Physics Branch of 
the Space Sciences laboratory has undertaken 
meteor observations using an image orthicon system. 
With this system, observations of larger mass 
meteoroids may be extended to as small as 10' 4 gm. 
The system that is used for meteor observations is a 
Maryland Telecommunications, Inc. Camera Chain 
with a GB 7967 low light level tube. The camera 
system uses a 525 line scan at the normal television 
frame rate of 30 frames per second. The optics 
used is an Aerojet Delft Rayxar lens with an aperture 
of approximately 5. 5 inches and a speed of f/0. 75. 

The results are video taped by Ampex 660 B/C tape 
recorders. 

Observational procedures are rather simple in 
that the entire system is portable and may be easily 
transported to different sites. Once in place, the 
camera is simply pointed to a selected region of the 
sky and the system is operated. Although our main 
observing site is presently on the south end of 
Redstone Arsenal, the system may be mounted in a 
DC-3 aircraft as shown in Figure 1 and flown to 
altitudes of 10 000 ft. This helps to alleviate the 
poor visibility which often prevails in the Huntsville 
area. With the Rayxar lens the field of view is 
approximately 13 by 16 degrees. Viewing at a 45- 
degree elevation, the observed area of the upper 
atmosphere in which meteor phenomena occur is 
approximately 1000 km 2 . 


Analysis of the observed results is a very 
painstaking effort. Many hours of analysis are 
required for each hour of observation, and a large 
backlog of data easily develops. The fact that pro- 
cedures are still being developed to more easily and 
accurately analyze the data also increases that 
backlog. 

The first step in the analysis is to examine each 
video tape thoroughly to locate all detectable meteors. 
Our present analysis procedures require that each 
observed area of sky be examined twice, since faint 
meteors are often difficult to detect. When a meteor 
is detected, its place on tape is noted along with 
information concerning its direction, spatial position, 
brightness, and speed. A meteor as observed by the 
image orthicon system is shown in Figure 2. 

Scanning the tapes results in meteor numbers 
and hourly rates. To date nearly 600 meteors have 
been detected. On the basis of the results recorded 
by other observers, hourly rates of from 15 to 20 
were expected. However, our analysis procedures 
have detected consistently higher rates than this; 
hourly rates often rise above 30 per hour and 
occasionally rise as high as 60 to 70 per hour. It 
should be mentioned that the higher rates were 
observed during late summer and fall of last year, 
and the highest rates were observed in the vicinity 
of the Leonid meteor shower. However, because of 
inclement weather , no observations were made during 
the actual shower period. 

The next step of the analysis is to determine the 
angular rates of the observed meteors. This con- 
sists of establishing two points along the trajectory 
of the meteor, finding the time required by the 
meteor to traverse this distance, and, with the aid 
of star charts, calculating the angular distance 
between the points. This results in the angular rate 
of the meteor in degrees per second. 

The distribution of angular rates as observed by 
the camera for 200 meteors is shown in Figure 3. 

This does not necessarily reflect the actual distribu- 
tion of angular rates, since it is easier for the 
orthicon to detect slower meteors. This results 
from the fact that a faster meteor has a faster writing 
speed across the photocathode, and fewer photons 
are deposited in each resolving element. Therefore, 
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Figure 1. Image orthicon and accessory equipment mounted in the 
DC -3 for airborne meteor observations. 


given two meteors of equal brightness, the slower 
object will appear brighter to the camera than the 
faster one. 

This effect is shown more clearly in Figure 4. 
In a laboratory test an object of steady brightness 
was observed at different angular velocities. It may 
be seen that the faster the object moved, the more 
difficult it was to detect. For meteors this means 
that as the angular rate becomes greater, the 
limiting magnitude becomes brighter. Limiting 
magnitudes for stars usually occur at about the 
9th or 10th magnitude, and in exceptional cases, at 
the 11th magnitude, but meteors this faint could only 
be detected if they possessed nearly a zero angular 
rate , and then the lack of a trail would probably 
render them unobserved. 


Determination of the magnitude of a meteor and 
the calculation of the integrated light curve of the 
meteor have presented the most difficult problems 
associated with the analysis to date. A complete 
analysis must take into consideration not only how 
to interpret what is actually observed, but also the 
manner in which the television system records what 
is observed. Low light level photometry through 
the use of television systems is a relatively new area 
of endeavor compared with the more frequent use of 
pattern discernment, and the analysis procedures 
are somewhat uncertain. A number of methods 
have been attempted to calculate the meteor magni- 
tudes, but no single method has shown itself able to 
determine these magnitudes over the range of 
observed brightnesses. It is possible that any 
complete analysis will have to rely upon more than 
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Figure 3. An angular rate distribution of meteors as recorded by the image orthicon 
during observations throughout November, 1961). 
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one method of brightness determination. Basically, 
the voltage output of a meteor will define its magni- 
tude when a comparison is made with stars of known 
magnitude panned across the field of view at known 
rates. 

At present, analysis is accomplished by viewing 
the television monitor with a photomultiplier tube. 

The face of the monitor is masked off except for a 
small slit through which the passage of the meteor 
may be observed. A peak-height light curve is then 
read from the photomultiplier tube by a storage 
oscilloscope. A typical light curve of a bright meteor 
is shown in Figure 5. Forty such light curves have 
been obtained. This method works for brighter me- 
teors, and meteor brightnesses have been determined 
to the seventh magnitude. Fainter meteors are lost 
in the background noise and must be examined by an 
alternate method. 

By the end of the summer, a second television 
system, an SEC vidicon, will complement our obser- 
vations. Using the two systems stationed at different 
sites, the meteor orbit, velocity, and height may be 
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Figure 5. A typical peak-height light curve as recorded by a 
photomultiplier tube observing the television monitor. 
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Figure 4. The decrease in apparent brightness of an 
object with increasing angular velocity. 
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calculated in addition to the information already investigations to be made into other areas of low 

recorded. The use of the two systems will also allow light level astronomy. 
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PHOTOELECTRIC OBSERVATIONS OE METEORS TO 6.4 MAGNITUDE 


By 

R. J. Naumann 


SUMMARY 

A photoelectric technique for detection and 
photometric measurement of meteors has been 
demonstrated. Through this technique meteors 
as faint as 6. 4 magnitude have been observed, 
although sky conditions at the present observing 
site limit the useful sensitivity to approximately 
5. 0 magnitude. The present observations indicate 
somewhat higher influx than expected by extrapola- 
tion of the photographic results, but the limited 
observations to date do not allow a determination 
of whether this result is caused by shower effects 
or by difference in calibration techniques. 


INTRODUCTION 

The definition of the meteoroid environment 
has been the object of a considerable effort within 
NASA. The presently accepted environment [1] 
depends heavily on photographic and visual observa- 
tions to define the influx of the larger masses. The 
observational limit of such techniques is approxi- 
mately fifth magnitude, but the useful limit, beyond 
which the probability of detecting a meteor is sub- 
stantially diminished, is approximately 2. 5 mag- 
nitude. This corresponds to approximately 0. 1 gm 
mass. The present NASA environment extrapolates 
from these data to the Pegasus data at 10“ 6 gm. 

The most critical region for large manned space- 
craft is in the mass range of 10“ 2 to 10“ 4 gm, which 
is exactly in the middle of the five orders of magni- 
tude extrapolation between existing data points. 

Since meteoroid impacts may represent the most 
serious hazard to long-term operation of large 
orbital space stations, it is imperative to refine 
our present knowledge of the meteoroid environment 
in this critical region. 

The most obvious approach is to extend the 
existing photographic data to fainter magnitudes by 
using more sensitive equipment such as extremely 
low light level TV systems. The Meteoroid Physics 


Branch of the Space Sciences Laboratory of Marshall 
Space Flight Center is currently using such a system 
that can observe meteors as faint as eighth magni- 
tude. However, such magnitudes can only be obtained 
for meteors with slow angular rates. This is because 
a fast meteor spends only a small fraction of the 
integration time of the instrument in a single 
resolving element. Consequently, the difference in 
output of that element is small compared with 
statistical fluctuations in surrounding elements, and 
the contrast is too low to allow detection unless the 
intensity is substantially above threshold. 

To supplement these TV observations and to 
aid in the photometry of meteors observed with 
video techniques, a 1. 5-m searchlight mirror is 
used to image a portion of the sky on the photo- 
cathode of an Amperex 150 UVP photomultiplier 
tube. The output of the photomultiplier tube is 
displayed on a memory oscilloscope through a type 
1A7A plug-in with the bandpass set for 0. 1 to 100 Hz. 
This was done to eliminate the dc background and the 
high-frequency shot noise from the sky. Meteors 
passing the field of view are seen as a momentary 
increase in signal that is stored on the memory 
scope and photographed by the observer. 

A similar technique was used by Davis [2] to 
determine color indices of meteors as faint as fourth 
magnitude. Davis’ system used a much larger 
120-mm photomuliplier tube with no optics. A 
50-degree field of view was defined by an aperture 
above the tube. Our system is more sensitive 
because of the smaller field and larger light- gathering 
capability. 


THEORETICAL LIMIT OF A PHOTOMULTIPLIER 
METEOR DETECTION SYSTEM 

The determining factor in the limiting detectability 
of the system is the fluctuations in the background. 
These fluctuations occur because energy is quantized, 
and the arrival of these quanta is a purely random 
process. In any random process, if N events are 
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expected in some interval of time t, the standard 
deviation is Wn (provided N » 1). In a photo- 
detector, photoelectrons are produced by photons 
within a certain energy range striking the photo- 
cathode. The number of electrons produced is 

where € is the quantum efficiency and is the 

number of photons within the energy range from the 
background. 

The average current is 
eN. e 

, b • _ 


If a randomly occurring signal is less than the 
noise, there is no way of detecting it. As the signal 
is made larger than the noise, the detection proba- 
bility is increased. A signal-to-noise rutio of five 
is generally considered sufficient to assure detection. 
This is equivalent to rejection of all but five sigma 
statistical fluctuations. 

Amplification will increase both signal and noise 
while adding additional noise. It will be assumed 
here that the added noise is negligible so that the 
system will be limited by the fundamental quantum 
noise at the input. This, of course, will set an 
upper limit on the performance of the system which 
no amount of design can surpass. 


* 

where e is the electronic charge and is the 
average number of photons per unit time. 

A standard deviation in i is 


(Ai) 


rms 


%/c We 
D 

t 


Let us now turn to the task of computing the noise 

background ft, and the signal ft . The rate at 
o s 

which signal photons strike the detector is 


N q = (collecting area) 


x (flux of photons within the proper 
spectral region) 


Since N, = N.t, 

D D 


(Ai) 


rms 



It is useful to describe the source in terms of stellar 
magnitude. By definition of stellar magnitude, the 
intensity of an Mth magnitude object is related to a 
zero magnitude object by 


Notice that the root mean square noise current can 
be reduced by increasing the measuring time t. 
This is equivalent to reducing the bandwidth of the 
system by excluding higher frequencies, since 
t= (2Af) _1 . 

Having defined the root mean square noise 
photocathode current, we now consider the signal 
current. If signal photons arrive at rate ft , the 
average signal current is 8 



The signal-to-noise ratio is 


I(m) = I(o) 10"°* 4M 

A zero magnitude star produces a spectral flux of 
3 x 10" 8 W/m 2 /\x on earth at a wavelength of 0. 5^t. 
Assume that this spectral flux is constant over a 
small AA centered about 0. 5/u. The energy flux 
from an Mth magnitude object within spectral range 
A A is 

I = 3 x 10“ 8 "°‘ 4M AX W/m 2 . 

To convert this to number of photons, we must 
divide by the average energy of a photon. A photon 
with a wavelength of 0. 5 p has an energy of 


SNR 


(Ai) 


sig 


rms 



_ JlL _ 6- 63 x 10" 34 x 3 x 10 8 
X 0. 5x iO- 6 

E = 4 x IO” 18 J 
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Therefore, 


Finally, 


N 

s 


n d 2 „ „„ Jrt ll-0.4M_ photons 

— 7 " x 0.75 x 10 AX “*■ 

4 sec 


N, 


7Td 2 

4 


(m 2 ) 


where AX is in microns, and d is the diameter 
of the objective lens in meters. 


X 0. 75 x 10 7 AX 


photons 


sec mr 10th mag star 


The background sky is an extended source and 
therefore the N. is 

D 


x N 


10 


10th mag star 
deg* i 


X (0.01745)^ degS 


= (collecting area) 
x (no. of photons/area/ solid angle) 
x (solid angle seen by detector) 


To extend to as faint objects as possible, it is 
evident that d must be as large as practical, e 
should be as large as possible, AX as large as 
possible, and N 10 as small as possible. Decreasing 
the solid angle will increase the sensitivity, and 
will also limit the number of events seen. This will 
be considered later. 


The steradiancy, or brightness, of the sky is often 
expressed in terms of a sky index, N 10 , the 
equivalent number of lOth-magnitude stars per square 
degree. From the result derived above, the number 
of photons in interval AX centered at 0. 5/x from a 
10th- magnitude star is 0. 75 x 10 7 AX photon/ sec/ m 2 . 


The present system uses a photomultiplier tube 
with broad response to increase AX and a large 
area searchlight mirror. Approximate values are 


d = 1. 5 m 


The solid angle subtended by the photocathode 
depends on the focal length and the photocathode 
area, 



F = 0.5m 


= _ir_ t0. 032) 2 
K 4 


nr 


€ = 0. 15 electrons/ photon 


By definition of solid angle and the fact that rays 
pass through the vertex of the lens undeviated, 



AX = 0. ipi 

N 10 = 2000 10th magnitude stars/deg 2 , 
The solid angle is 


Since 1 degree is 0.01745 radian, one square degree 
subtends (0. 01745) 2 . 


7r (0. 032) 2 
4. 5* (0. 01745) 2 


10. 5 deg 2 
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The area of the collector is 

Area = = 1.76 m 2 

4 


The ultimate performance of this technique may 
be estimated by considering possible improvements. 
The expression for limiting magnitude is 


N, = 1. 76 x 0. 75 x 10 7 x 0. 1 x 2000 x 10. 5 
b 


M 


Lim 


const + 1. 25 log 


€d 2 A\ 

N 10 nAf 


= 2. 8 x 10 10 photons/ sec 

• 11-0. 4M 

N = 1.76 x 0.75x 10 x 0. 1 

s 

j on 10-0. 4M . 

- 1.32x10 photons/sec 

i = 2.8x 10 19 x 15 x 1. 6 X 10- 19 
av 

= 0. 67 x 10.- 9 A 


Not much improvement can be expected in e or 
AX. The collecting area is as large as can be 
obtained for a reasonable cost. The biggest 
potential gain lies in the sky background, N 10 . 

The value of 2000 was used for the vicinity of 
Huntsville. An improvement factor of four could 
be realized by moving the system to a darker site 
which would increase the limiting magnitude by 
0. 75. The solid angle could be reduced by using a 
smaller photocathode or a longer focal length. 
However, this is accompanied by a reduction in 
data acquisition rate since a smaller area is viewed 
and by increased bandwidth requirements if the 
meteor duration is limited by the field of view. 


(Ai) = 1. 6 x 10“ 19 -Vo. 15 x 2. 8 x 10 10 2Af 
rms 


The transit time of a meteor across the field 
of view is 


= 1.47 x 10 -1< n/HTa 


With a bandpass of 100 Hz, the root mean square 
background is 1. 47 x 10“ 13 A at the photocathode. 
The faintest meteor that could be detected against 
this background is obtained by requiring 


SNR = N /— — = 5 

8 J 2AfN 

b 

or 


I.MXIO 10 - 0 ' 4 " ■ 5 2x1.0 

V 0. 15 



where d^ is the diameter of the photocathode, F 

is the focal length, c o is the angular rate, and r 
is the smallest distance from the center of the 
photocathode to the image of the meteor trail. 



IMAGE OF 
METEOR TRAIL 


1.32 x 10 


10-0. 4M 


30. 6 x 10® 


Taking logarithms, 

10-0.4M = log 23.2 + 6 = 7.364 



0.4M = 2.636 


M 


Lim 


6.6 


where Ar is the area included by radius r. Setting 
this radical to a value of 0. 1 assures that 99 percent 
of all meteors within the field will have transit 
times greater than 
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t = 


d, 0. 1 
k 

Fa) 


3.2 

500u) 


6.4 

it) 


CALIBRATION 


For o) =0.7 rad/sec, the shortest time that 
max 

need be considered is 9 ms. The minimum bandwidth 
is ~ 50 Hz. 


Since the system is always operated in a region 
where the anode current is small enough to assure 
linearity, the voltage observed can be related to 
the magnitude of a meteor by 


The product S7Af must therefore be considered 
together. Since 


V- V L 


= (V -V.)10 
o b 


-0.4M 


Af = 


and 


1 

2 t 


1 


F 0.7 

2d, 0.1 
k 


3. 5 F 


where V - is the observed signal above back- 


ground, V 


is the voltage associated with a 


zero magnitude meteor above background, and M 
is the magnitude of the observed source. The 
quantity V q - is obtained by calibration against 

known stars. 


then 


£2Af 


^k 

F 


Thus sensitivity can be improved by 0. 375 magnitudes 
by reducing c^/F by two, but the field of view and 

data acquisition rate are reduced by a factor of four. 

The practical limit would then appear to be the 
present system using an il-mm photomultiplier, a 
bandwidth of 100 Hz 1 (which will still accommodate 
97 percent of even the fastest meteors) in a location 
with sky index of 500. This will gain 

0. 75 magnitudes from better sky 
1.16 magnitude from smaller solid angle 


1. 91 total gain 
6.60 existing system 


M 


Lim 


= 8.51 


In practice, since the background sky in the 
3-degree field of view is approximately the same 
as a zero magnitude star, the calibration star 
should be as bright as possible to reduce errors in 
obtaining the difference V - V^. Sirius is an 

ideal calibration star because of its -1. 58 magnitude 
and also because it is a type A1 star for which 
photographic and visual magnitude nearly coincide. 
For this reason, no correction will be made for 
color index, and all meteor magnitudes will be 
referred to type Al spectral class. 

A calibration for a typical observing night is 
described. A recording, integrating digital 
voltmeter was used to minimize fluctuations from 
twinkling. An average of 19. 364 V was obtained 
from Sirius. As Sirius moved out of the field, the 
reading dropped to 3. 890 V. Thus, 


19.364 - 3.890 


(V " V 10 
o b 


-0. 4 (-1. 58) 


V - V. = 
o b 


3.60 


The area will be reduced by a factor of nine, but 
this will be approximately compensated for by the 
increase in meteor population for decreasing size. 


Several other stars were measured to check this 
calibration. The results are summarized in the 
following table: 


1. For very faint meteors the duration may be limited by the meteor being consumed rather than by the field 
of view. For this reason 100 Hz is the smallest bandwidth that should be considered. 
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Star 

Type 

M 

V 

V . 
obs 

v b 

M (calculated) 

Capella 

G8 

0,21 

6.40 

3.51 

0. 24 

Procyon 

F5 

0.48 

5.92 

4. 12 

0.753 

Betelgeuse 

M2 

0.92 

5. 06 

4. 12 

1.58 

Pollux 

K0 

1.21 

4.47 

4. 03 

2.28 

Castor 

A1 

1.58 

1 

4.69 

4. 03 

! 1.84 


The M (calculated) was found from 


M cal = ** 39 " 2 * 5 log (V " V b J 


There are several possible reasons for the 
differences between M and Probably the 

major source of error lies in obtaining the sky 
background. For these stars the mirror was 
moved to a region adjacent to the star, and the 
voltage was read. A slight difference in the sky 
background will make a substantial error for the 
fainter magnitudes. Also, the different spectral 
classes can be expected to require increasing 
correction for increasing redness since the tube 
used has an S-13 UV response. 

The sky index can be obtained by setting 

V - V. = V * The V. , , . ranged from 

b b(obs) b(obs) & 

3. 2 V at zenith to 4. 12 V at an elevation of approxi- 
mately 30 degrees. This corresponds to 0. 128 
magnitude at zenith and -0. 145 magnitude at 
30-degrees elevation. In terms of 10th- magnitude 
stars/deg 2 , assuming a field of 10. 5 deg 2 , the sky 
index is 850 at zenith and 1100 at 30 degrees. 

These are probably low for the sky in the vicinity of 
Huntsville, where fourth-magnitude stars are barely 
visible to the eye. The response of the system falls 
off with angle before the geometric limit is reached 
because of spherical aberration. The exact shape 
of the response curve has not yet been determined, 
but from observing stellar transits, the response 
appears to be parabolic in nature with the vertex on 
the optical axis and the intercept at the geometric 
limit. It is estimated that the effective solid angle 
is 50 to 67 percent of the geometric solid angle. 



IDEAL GEOMET- 
RIC RESPONSE 
ESTIMATED AC- 
TUAL RESPONSE 


This correction will increase the sky index to 1500 
through 2000. 


A substantial portion of the sky background is 
backscatter from Huntsville lights. This is evident 
from the 120-Hz component in the sky background 
that is typical of mercury lamps (Fig. 1). This is 
averaged out in the calibration procedure by using 
an integration time of 1/60 sec, but, as will be 
seen later, this component severely hampers the 
detection of faint meteors. 


METEOR OBSERVATIONS 


For the preliminary results described in this 
report, meteors were detected by an observer 
monitoring the memory scope on slow scan. Typical 
meteors are shown in Figure 2. The broad band of 
noise of the 120-Hz component is caused by 
Huntsville mercury vapor lamps. Brighter meteors 
show up as a temporary level shift. However, 
fainter meteors of short duration can be detected 
only if their peak intensity coincides with a peak in 
the background. This caused many of the fainter 
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Figure 1. 120-Hz sky background and level 

shift caused by a 3.45 magnitude meteor 
(50 mV/cm, lOms/cm). 



Figure 2b. 5. 16-magnitude meteor 

(50 mV/cm, 50 ms/cm). 

meteors to be missed, as can be soen in the 
histogram (Fig. 3). Apparently, detection probabil- 
ity diminishes beyond fifth magnitude. This, of 
course, can be rectified by observing from a site 
with better viewing conditions and by using electronic 
techniques to trigger a counter when an event occurs 
that is 5X background noise. The goal for the 
present phase was to demonstrate that the technique 
could be used to extend optical meteor data to 
fainter magnitudes. 

The results of 18 hours of observation are listed 
in Table 1. The voltages associated with the meteors 
have been corrected for filter attenuation (Fig. 4), 
which was obtained experimentally by feeding pulses 
of various shapes and duration into the system and 
observing the ratio of the observed pulse to input 
pulse. The voltages were then converted to magni- 
tude using the calibration result 


M = 1. 39 - 2. 5 log 10 (V - V ) 



Figure 2a. Typical meteors observed — 
3. 58- magnitude meteors 
(50 mV/cm, 100 ms/cm). 



Figure 2c. 4. 5-inagnitude meteor 
(50 mV/cm, 50 ms/cm) . 

The results of these observations are plotted 
along with those of Hawkins and Upton [3J in 
Figure 5. The area was obtained by taking the 
portion of sky included in 10. 5 deg 2 at a distance of 
100 km. 

It is interesting to note that the present study 
compares reasonably with the extrapolated results 
of Hawkins and Upton. The slope indicated by the 
observed data is practically identical with that 
deduced by Hawkins and Upton. The present obser- 
vations indicate an influx rate about a factor of three 
to four higher than expected from extrapolating the 
Hawkins and Upton curve. It is not yet certain that 
this is real. The present observations were made 
over a limited time (approximately 1 month), and 
it is known that fluctuations occur from month to 
month. The observations also contain the time of 
the Geminids shower which may have added signifi- 
cantly to the number observed. At present, limited 
statistics do not allow a determination of the effect 
of showers. There is also the problem of definition 
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Figure 3. Results of this study compared with photographic work of Hawkins and Upton. 



Figure 4. Filter attenuation at different pulse durations. 
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TABLE 1. METEOUS OBSERVED IN 17 HOURS, 22 MINUTES OBSERVING TIME 


No. 

Time 

dc Level 
V 

Pulse Height 
mV 

Duration 

msec 

Magnitude 

H 

11-20-68 

12:45 

4.5 

187 

150 

3.22 

2 

11-21-68 

1:25 

4. 15 

150 

100 

3.45 

3 

1:34 

4. 00 

90 

94 

4.01 

4 

2:12 

3. 85 

41 

80 

4. 86 

5 

2:45 

3. 85 

12 

34 

4. 86 

6 

11-22-68 

12:05 

3.0 

1’ 

500 

150 

2. 14 

7 

12:36 

4. 1 

160 

80 

3. 58 

8 

12:40 

4. 1 

170 

GO 

3.32 

9 

12:45 

4. 1 

57 

380 

4. 50 

9a 

12:45 

4. 1 

36 

26 

5.00 

10 

1:02 


75 

52 

4. 20 

11 

1:50 

3.0 

13 

13 

6. 11 

12 

1:58 

3. 0 

322 

80 

2. 62 

13 

12-13-68 

10:05 

3. 7 

95 

42 

3. 95 

14 

11:10 

4. 6 

25 

226 

5.40 

15 

10:40 

3.7 

171 

19 

3.27 

16 

10:45 

3. 7 

40 

226 

4. 89 

17 

10:47 

3.6 

10 

42 

6.39 

18 

11:03 

3. 5 

70 

17 

4. 28 

19 

11:15 

3. 5 

208 

47 

3.09 

20 

12:00 

4.2 

62 

80 

4.41 

21 

12:15 

3. 9 

60 

316 

4.44 

22 

12:32 

4. 1 

82 

190 

4. 11 

23 

12:45 

3. 9 

102 

270 

3. 87 

24 

12-16-68 

9:00 

5.2 

107 

26 

3. 82 

25 

9:03 

5. 2 | 

35 

231 

5.03 

26 

9:05 

5. 0 

60 

35 

4.44 

27 

9:36 

4.4 

43 

10 

4. 81 

28 

9:37 

4.4 

40 

210 

4. 89 

29 

9:40 

4. 4 

71 

42 

4.26 

30 

9:41 

4.4 

70 

37 

4.28 

31 

12:20 

4.7 

50 

90 

4.64 

32 

1-13-69 
(9:05 - 12:00) 
8:56 

4. 01 

27 

80 

5.31 

33 

10:03 

3. 87 

31 

150 

5. 16 

34 

11:06 

3.62 1 

41 

21 

4. 86 
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Figure 5. Influx rates of this study compared with those observed by Hawkins and Upton. 


of magnitude since our system has a response 
different from that of photographic film. 

CONCLUSIONS 

The use of a photomultiplier tube and a large 
objective mirror appears to be a useful technique 


for counting and obtaining photometric data for 
faint meteors. The present study extended the 
Hawkins and Upton results from 2.5 to 5.0 
magnitude. Meteors near the theoretical limit 
of 6.6 were observed but with low detection 
probability. The technique appears to have an 
ultimate capability of 8.51, which corresponds to 
masses of 4 x 10 3 4 gm. 
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OBSERVATORY PLANS AND ASTRONOMICAL APPLICATIONS 

By 

James McGuire 


The observatory constructed at the southern end 
of Redstone Arsenal was developed to improve the 
meteor observation program at MSFC and to aid in 
developing other astronomical programs. The area 
for the site, shown in Figure 1, was selected be- 
cause of its favorable location and the good visibility 
conditions present. Table 1 shows that the observa- 
tional conditions in the Huntsville area are more 
desirable than at 50 percent of the other observatory 
sites in the United States. The data in this table 
were obtained from the U. S. Weather Bureau and 
a Brown Engineering Company report. 



Figure i. Location of observatory site at 
Redstone Arsenal. 


The views of the horizon observed from the site 
are shown in Figures 2 through 5. Figure 2 is a view 
of the eastern horizon; Figure 3, the southern hori- 
zon; Figure 4, the western horizon; and Figure 5 is 
a view of the northern horizon. The building shown 
in Figure 2 is not a part of the observatory, but 
belongs to the Army. The site is atop a hill approxi- 
mately 300 feet tall, the same height as the vibration 
test tower pictured in the background of Figure 5. 


TABLE 1. OBSERVATIONAL CONDITIONS 
OF OBSERVATORY SITE 


Sky Glow: 

40 percent better than 
Monte Sano; 40 percent 
worse than the best loca- 
tion within 30 miles of 
Huntsville. 

Resolution Limitation: 

0. 15 second of arc 

Clear Nights: 

15 percent better here than 
for the average eastern and 
midwestern observatory; 

15 percent worse here than 
for the average western 
observatory. 


Two telescopes and two domes will be used when 
the site becomes operational. The domes are made 
by Ash-Domes Corporation; the larger is 18 feet in 
diameter and the smaller, shown in Figure 6, is 12 
feet in diameter. One of the telescopes is manufac- 
tured by Ferson and has a 16- inch primary mirror 
(Fig. 7). This model can be used as a Newtonian 
telescope with a 53-minute field of view or as a 
Cassegrain telescope with a 15-minute field of view. 
The second telescope, shown in Figure 8, was made 
by Celestron Pacific and is also a 16 inch. This 
telescope is similar to a Schmidt telescope since it 
has a correcting lens for the spherical primary 
mirror. The advantage of this type of telescope is 
that it is a closed system and is therefore easier to 
keep clean and is not affected by heat waves in the 
telescope tube. When used in the Cassegrain mode, 
this telescope has a 20-minute field of view. The 
Celestron telescope also has a rich field adaptor 
that will allow a 1-degree field of view. 

The instruments to be used on the telescope will 
be low light level TV systems, photographic equip- 
ment in conjunction with image intensifiers, photo- 
metric equipment, and, possibly, a spectrograph. 
Figure 9 shows a Celestron telescope with a low 
light level TV camera attached and the rich field 
adaptor being used. 


Ill 
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UNITS OF MEASURE 


In a prepared statement presented oi* August 5, 1965, to the 
U. S. House of Representatives Science and Astronautics Committee 
(chaired by George P, Miller of California), the position of the 
National Aeronautics and Space Administrationon Units of Measure 
was stated by Dr. Alfred J. Eggers, Deputy Associate Administrator, 
Office of Advanced Research and Technology: 

"In January of this year NASA directed that the international 
system of units should be considered the preferred system of units, 
and should be employed by the research centers as the primary 
system in all reports and publications of a technical nature, except 
where such use would reduce the usefulness of the report to the 
primary recipients. During the conversion period the use of cus- 
tomary units in parentheses following the SI units is permissible, 
but the parenthetical usage of conventional units will be discontinued 
as soon as it is judged that the normal users of the reports would 
not be particularly inconvenienced by the exclusive use of SI units. " 

The International System of Units (SI Units) has been adopted 
by the U. S. National Bureau of Standards (see NBS Technical News 
Bulletin, Vol. 48, No. 4, April 1964). 

The International System of Units is defined in NASA SP-7012, 
"The International System of Units, Physical Constants, and 
Conversion Factors,” which is available from the U. S. Government 
Printing Office, Washington, D. C. 20402. 

SI Units are used preferentially in this series of research re- 
ports in accordance with NASA policy and following the practice of 
the National Bureau of Standards. 
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